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Abstract 

Tempest  is  a  computational  tool  that  is  being  developed  to  predict  the  seakeeping  and 
dynamic  stability  performance  of  a  steered  ship  in  large  waves.  The  theory  for  Tempest  is  being 
developed  in  several  “levels”.  The  difference  between  the  theory  levels  is  in  the  fidelity  and 
complexity  of  the  environment  definition  and  the  component  force  models.  This  report  describes 
the  Level-0  theory.  The  objective  of  the  Level-0  version  of  Tempest  is  to  establish  a  working 
computational  tool  that  provides  reasonable  seakeeping  predictions  and  to  establish  the 
framework  of  the  Tempest  tool  into  which  the  more  accurate  Level-II  and  Level-Ill  models  can 
be  implemented.  The  Level-0  theory  is  based  on  simple  models  that  include  a  linear 
hydrodynamic  force  model  based  on  user-supplied  added  mass  and  damping  coefficients,  a 
nonlinear  model  for  the  hydrostatic  and  Froude-Krylov  forces,  and  empirical  models  for  the 
maneuvering  forces  and  the  forces  from  the  bilge  keels,  rudders  and  propeller.  This  report  is  a 
compilation  of  the  equations  of  motion,  wave  environment,  and  component-force  theory  white 
papers  provided  to  the  Tempest  code  developers,  DRS  Defense  Solutions,  for  implementation  of 
the  Level-0  theory. 


Administrative  Information 

The  work  described  in  this  report  was  perfonned  at  the  Naval  Surface  Warfare  Center, 
Carderock  Division  (NSWCCD)  by  the  Seakeeping  Division,  Code  5500.  The  development  of 
the  Tempest  Level-0  theory  white  papers  described  within  this  report  was  funded  by  Naval  Sea 
Systems  Command  (NAVSEA),  PMS500,  under  Program  Element  0604300N.  The  work  was 
performed  predominately  under  work  unit  07-1-5500-752-10.  The  writing  of  this  report  was 
funded  by  NAVSEA  05D1  under  the  CPSD  SE/TA  program  supporting  Tempest  development. 
The  associated  Program  Element  is  0603563N  and  work  unit  is  09-1-2124-206. 


Introduction 

Tempest  is  a  computationally  efficient  time  domain  tool  designed  to  predict  the  motions  of 
a  steered  ship  in  large  waves.  It  is  intended  for  use  as  both  a  seakeeping  and  a  dynamic  stability 
prediction  tool.  The  planning  for  the  development  of  Tempest  was  initiated  in  2006.  The  plan 
called  for  a  multi-phased  effort  that  would  promote  early  progress  on  the  code  development 
while  the  theory  is  in  development,  with  the  goal  being  that  a  more  polished  and  complete 
program  would  be  available  sooner.  The  first  phase  of  the  code  development  is  focused  on 
creating  the  foundational  code  across  all  three  program  tiers  (Client,  Business,  and  Data)  and  is 
therefore  not  expected  to  be  used  for  actual  simulation  work  applied  to  real  ship  designs. 
However,  while  there  is  not  a  need  to  have  accurate  simulation  results,  the  components  of  the 
ship  motion  theory  need  to  be  included  to  the  greatest  extent  possible  in  order  to  provide  a 
meaningful  basis  for  developing  code  related  to  the  inputs  and  outputs  of  the  Client  tier,  the 
interactions  of  the  integrators,  force  models,  and  other  foundational  aspects  in  the  Business  tier, 
and  finally  the  persistence  capabilities  and  needs  in  the  Data  tier.  To  support  this  need,  a  so- 
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called  “Level-0”  theory  has  been  established  that  provides  a  nearly  complete,  though  lower 
fidelity,  ship  motion  theory  to  be  implemented  in  the  first  release  of  Tempest.  The  Level-0  code 
is  mainly  based  on  simple  theories  that  can  be  used  to  develop  a  working  code  that  produces 
“reasonable”  predictions. 

Several  “levels”  of  theory  will  be  developed  for  Tempest.  All  levels  of  the  theory  are 
based  on  the  same  rigid  body  equations  of  motion  and  time  stepping  technique.  The  difference 
between  the  levels  is  in  the  fidelity  and  complexity  of  the  environment  definition  and  the 
component  force  models.  This  theory  manual  describes  the  theory  implemented  in  “Level-0” 
Tempest  code.  The  objective  of  the  “Level-0”  theory  is  to  develop  a  well-documented 
framework  into  which  the  more  accurate  Level-II  and  Level-Ill  theories  can  later  be 
implemented.  The  Level-0  theory  will  use  very  approximate  methods  for  certain  force 
components  such  as  the  radiation,  diffraction  and  maneuvering  force  models,  as  these  will  be 
replaced  by  completely  new  force  models  as  the  Level-II  and  Level-Ill  theories  are  developed. 
Certain  components  of  the  Level-0  theory  such  as  the  equations  of  motion  and  time  integration 
schemes  and  the  hydrostatic  and  Froude-Krylov  force  models  are  intended  to  carry  over  to 
Level-II  and  possibly  even  into  Level-Ill  with  only  minor  modifications. 

The  theory  described  in  this  document  was  developed  from  a  variety  of  sources.  The 
outline  of  the  theoretical  structure  of  Tempest  was  developed  by  a  core  group  at  Naval  Surface 
Warfare  Center,  Carderock  Division  (NSWCCD).  A  Theory  Advisory  Panel  (TAP)  was  formed 
to  review  and  assist  in  the  theory  development  for  Tempest.  A  series  of  white  papers  were 
written,  with  each  white  paper  discussing  a  specific  aspect  of  the  code. 

The  Level-0  theory  has  been  developed  in  the  Seakeeping  Division  (Code  55)  of  NSWCCD 
in  West  Bethesda,  MD.  Code  development  has  been  performed  by  DRS  Defense  Solutions  in 
Stevensville,  MD.  The  theory  has  been  documented  and  communicated  to  the  code  developers 
incrementally  in  the  fonn  of  “Level-0  theory  white  papers”. 

The  present  document  provides  a  record  of  these  theory  white  papers  within  a  single  report 
that  is  able  to  be  referenced.  The  following  document  is  not  intended  to  be  a  comprehensive 
theory  manual,  but  rather  a  compilation  of  the  theory  documents  from  which  the  Level-0  code  is 
developed. 
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1  Nomenclature 

The  nomenclature  tables  below  describe  the  variables  used  in  the  Tempest  theory 
development.  As  this  is  a  theory  manual  that  has  been  developed  from  the  compilation  of 
several  component  theory  documents,  the  nomenclature  has  been  broken  into  several  sections 
corresponding  to  the  significant  theory  areas.  The  first  section  describes  the  nomenclature 
pertaining  to  the  wind  and  wave  environment.  The  second  section  describes  the  nomenclature 
used  for  the  general  equations  of  motions  and  the  coordinate  systems  used  throughout  the  theory. 
The  remaining  sections  describe  the  nomenclature  used  in  the  various  force  modules. 

The  nomenclature  will  be  consistent  throughout  the  theory  manual.  In  some  instances  a 
symbol  may  have  multiple  definitions.  For  instance  the  symbol  “7”  may  mean  wave  period  in 
the  wave  and  wind  environment  theory,  thrust  in  the  propeller  force  module  theory,  and  draft  in 
the  bilge  keel  damping  theory.  In  these  instances  it  will  always  be  obvious  which  meaning  the 
symbol  has  in  any  given  equation.  Where  possible  these  conflicts  are  avoided  by  not  using  the 
same  symbol  for  different  variables,  but  only  to  the  extent  to  which  it  helps  to  clarify  the  theory. 
If  a  commonly  accepted  symbol  exists  (i.e.,  T  for  both  period  and  draft,  X  for  wavelength,  etc.) 
then  that  symbol  is  always  used. 

Careful  attention  is  paid,  however,  to  make  sure  that  two  different  symbols  are  not  used  for 
the  same  variable  (i.e.  using  both  Q  and  77  for  wave  elevation  in  different  sections  of  the  theory 
manual).  With  regard  to  the  coordinate  system  definitions,  upper  case  ( X,  Y,Z)  are  used  to  refer  to 
coordinates  in  the  earth-fixed  system  and  lower  case  (x,y,z)  are  used  to  refer  to  coordinates  in  the 
ship-fixed  system. 


1.1  Wave  Environment 


Svmbol 

Definition 

SI  Units 

«/' 

(Real-valued)  first  order  wave  amplitude  equal  to  half  the 
distance  from  the  trough  to  the  crest 

m 

Ojk 

(Real-valued)  second  order  wave  amplitude 

m 

c 

Wave  celerity  (phase  velocity) 

m/s 

cg 

Group  velocity  of  wave 

m/s 

E() 

Expected  value  of  a  random  variable 

f 

Wave  frequency  in  units  of  cycles  /  second 

cycles/s 

g 

Acceleration  due  to  gravity 

m/s2 

G 

Directional  spreading  function 

h 

Water  depth 

m 

Hs 

Significant  wave  height 

m 
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Wave  number 


rad/m 


nth  moment  of  a  spectrum 
Number  of  wave  components 
Number  of  spreading  angles 
Spectral  density 
Time 

Wave  period 
Modal  wave  period 

Wave  orbital  velocities  in  earth- fixed  X,  Y  and  Z  directions 

Wave  orbital  accelerations  in  earth- fixed  X,  Y  and  Z 
directions 

Modified  Z  position  for  pressure  stretching 
Quadratic  transfer  function 
Wheeler’s  depth  decay  modifier 

Direction  of  wave  propagation  measured  from  the  positive 
X-axis  counterclockwise  about  the  Z-axis 

Direction  of  primary  wave  system 
Spreading  angle 

Lower  limit  on  wave  component  direction  for  spread  seas 

Upper  limit  on  wave  component  direction  for  spread  seas 

Direction  of  wind  measured  from  the  positive  X-axis 
counterclockwise  about  the  Z-axis.  If  ship  is  traveling  in  the 
positive  X  direction,  0  will  be  a  tail  wind. 

Wave  phase  angle 
Velocity  potential 

JONS  WAP  peak  enhancement  factor 
Wave  length 
Water  density 
Wave  frequency 

Lower  truncation  limit  of  the  spectral  frequencies 
Modal  wave  frequency 


m2(l/sec)n 

nrsec/rad 

s 

s 

s 

m/s 

m/s2 

m 

radians 

radians 

radians 

radians 

radians 

radians 

radians 

m2/s 

m 

kg/m3 

rad/s 

rad/s 

rad/s 


Upper  truncation  limit  of  the  spectral  frequencies 

rad/s 

C 

Wave  elevation 

m 

1.2  Equations  of  Motion  and  Coordinate  Systems 

Symbol 

Definition 

SI  Units 

Aijoo 

Added  mass  at  infinite  frequency  in  the  zth  mode  due  to  unit 
acceleration  in  the  jth  direction 

kg 

V 

Added  mass  at  wave  encounter  frequency 

kg 

F 

Total  force  vector  in  ship-fixed  reference  frame 

N 

F' 

Total  force  vector  in  ship-fixed  reference  frame,  not 
including  weight  and  infinite  frequency  added  mass  force 

N 

1} pc?  lyyi  Izz 

Moments  of  inertia  about  ship-fixed  origin 

kg  nr 

Fyt  lyzi  Izx 

Products  of  inertia  about  ship-fixed  origin 

kg  nr 

M 

Total  moment  vector  in  ship-fixed  reference  frame 

Nm 

M ' 

Total  moment  vector  in  ship-fixed  reference  frame,  not 
including  moment  due  to  ships  weight  and  infinite  frequency 
added  mass 

Nm 

m 

Mass  of  the  ship 

kg 

P,  q,  r 

Angular  velocities  about  the  three  ship-fixed  axes 

rad/s 

p,q,r 

Angular  accelerations  about  the  three  ship-fixed  axes 

rad/s' 

U,  V,  w 

Velocities  in  the  x,v,z  directions  (ship-fixed),  respectively 

m/s 

u,v,w 

Accelerations  in  the  x,y,z  directions  (ship-fixed). 

m/s2 

W 

Weight  of  the  ship  (m-g) 

N 

X,Y,Z 

Position  of  the  origin  of  the  ship-fixed  axes  system  relative  to 
the  earth-fixed  axes  system 

m 

x,y,z 

Velocity  of  ship  origin  in  the  earth- fixed  frame 

m/s 

&1 1  &2n  & 3 9  & 4 

The  components  of  the  quaternion 

cj) ,  9 ,  y/ 

Euler  angles:  roll,  pitch  and  yaw 

rad 

<f> ,  9 ,  y/ 

Time  derivatives  of  Euler  angles 
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1.3  Bilge  keel  force  module 


Symbol  Definition  SI  Units 

A(/t  Roll  amplitude  rad 

bBK(x)  Bilge  keel  breadth  (height)  m 

C„  Block  coefficient 

D 

CM  Midship  section  coefficient 

Lbk  Bilge  keel  length  m 

L  Length  between  perpendiculars  m 

OG  Vertical  center  of  gravity  w.r.t.  waterline  m 

P  Roll  velocity  rad/s 

rBK  (x)  Mean  distance  from  G  to  bilge  keel  m 

R  Bilge  radius  m 

T ( x )  Sectional  draft  m 

U  Ship  speed  m/s 

Ul0  Ship  reference  speed  =  lOm/s  m/s 

VCG  Vertical  center  of  gravity  w.r.t.  keel  m 

x  Longitudinal  coordinate  m 

5  Angle  between  the  line  connecting  CG  and  the  bilge  keel  rad 

root  and  the  horizontal  through  CG 

Y  Angle  between  the  line  connecting  CG  and  the  bilge  keel  rad 

root  and  the  bilge  keel  plane 

P  Density  of  water  kg/m 

cr(x)  Sectional  area  coefficient 

co^  Roll  frequency  rad/s 
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1.4  Rudder  Force  Module 


Symbol 

Definition 

SI  Units 

a 

Aspect  ratio  of  rudder  (geometric  aspect  ratio) 

ae 

Effective  aspect  ratio 

Ar 

Area  of  rudder  planform 

nr 

Cr 

Chord  at  root  of  rudder 

m 

Ct 

Chord  at  tip  of  rudder 

m 

Xc>yc>Zc 

Location  of  the  center  of  pressure  of  the  rudder 

m 

Xm,ym,Zm 

Location  of  the  midspan  point  along  quarter  chord  line 

m 

Xr<)'  r<- r 

Location  of  the  rudder  root  quarter-chord  point 

m 

xt.yt.zt 

Location  of  the  rudder  tip  quarter-chord  point 

m 

b 

Mean  rudder  span.  Note  that  this  is  sometimes  referred  to  as 
the  “semi-span”  such  as  in  Whicker  and  Lehlner,  1958 

m 

C~ 

Mean  rudder  chord  (average  of  tip  and  root  chord) 

m 

CL 

Lift  Coefficient 

Cd 

Drag  Coefficient 

Cdc 

Cross  flow  drag  coefficient 

Cd0 

Minimum  section  drag  coefficient  (0.0065  for  NACA  0015 
sect.) 

CN 

Normal  force  coefficient  for  a  wing  in  stalled  condition  as 
described  in  Hoerner  (1975).  Recommended  value  1.8  to 

2.0. 

FR 

Lorce  vector  for  rudder  force  in  ship-fixed  reference  frame 

N 

L 

Lift  from  rudder 

N 

mr 

Moment  vector  defining  moment  from  rudder  in  ship-fixed 
reference  frame. 

Nm 

d 

Drag  on  rudder 

N 

D 

Diameter  of  propeller 

m 

T 

Propeller  thrust 

N 

P,  q,  r 

Angular  velocities  about  the  three  ship-fixed  axes 

rad/s 

u,  v,w 

Velocity  of  the  ship  at  the  ship-fixed  reference  frame  origin 
in  the  x,y,z  directions  (ship-fixed),  respectively. 

m/s 

UWr)  Vw-, 

Wave  orbital  velocity  at  the  rudder,  defined  in  the  ship-fixed 

m/s 
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frame 

wf  Wake  fraction  coefficient 

a  Angle  of  attack  rad 

as  Stall  angle  rad 

fj  Drift  angle  of  flow  into  rudder  rad 

8  Rudder  deflection  angle  rad 

3 

p  Density  of  water  in  which  rudder  is  operating  kg/m 

X  Taper  ratio  (ct/cr) 

A  Sweep  angle  of  quarter-chord  line  rad 


yR  Flow  straightening  factor,  used  to  adjust  the  angle  of  attack 

of  the  rudder  to  account  for  the  flow  straightening  influence 
of  the  hull  and  propeller,  (default  value  1 .0) 

(f>R  Dihedral  angle  of  rudder,  which  is  measured  as  the  rad 

inclination  from  vertical  with  a  positive  angle  indicating  the 
rudder  tip  is  further  to  starboard  than  the  rudder  root. 


1.5  Propeller  Force  Module 


Svmbol 

Definition 

SI  Units 

Ad 

Area  of  propeller  disk 

m2 

As 

Area  of  submerged  portion  of  propeller  disk 

m2 

D 

Diameter  of  the  propeller 

m 

F  Cross 

Component  of  the  propeller  force  perpendicular  to  the 
propeller  shaft  in  the  direction  of  the  cross  flow 

N 

F Norm 

Component  of  the  propeller  force  perpendicular  to  both  the 
propeller  shaft  and  cross  flow 

N 

FP 

Force  vector  for  propeller  force  in  ship-fixed  reference  frame 

N 

J 

Advance  coefficient 

- 

Kj 

Propeller  thrust  coefficient 

- 

MP 

Moment  vector  defining  moment  from  propeller  in  ship-fixed 
reference  frame 

Nm 
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n  Rotational  speed  of  the  propeller,  revolutions  per  second  rev/s 

ns  Unit  vector  defining  direction  of  propeller  shaft  in  ship-fixed 

frame 

nc  Unit  vector  defining  direction  of  cross  flow 

iiN  Unit  vector  defining  the  direction  normal  to  cross  flow  and 

propeller  shaft 

p,  q,  r  Angular  velocities  about  the  three  ship-fixed  axes  rad/s 

R  Propeller  radius  m 

RHP  Integer  flag  defining  whether  propeller  is  right  hand  turning 

or  left  hand  turning  (+1  for  right  hand  turning,  -1  for  left 
hand  turning) 

T  Propeller  thrust  (component  of  propeller  force  parallel  to  N 

propeller  shaft) 

tp  Thrust  deduction  coefficient 

Ua  Component  of  velocity  vector  at  center  of  propeller  disk  that 


is  parallel  to  the  propeller  shaft 

Up  Total  velocity  vector  at  the  center  of  the  propeller  disk  m/s 

77  Cross  flow  velocity  vector  at  the  center  of  the  propeller  disk  m/s 

V A  Advance  velocity  into  propeller  disk  m/s 

u,  v,  w  Velocity  of  the  ship  at  the  ship-fixed  reference  frame  origin  m/s 


in  the  x,y,z  directions  (ship-fixed),  respectively 

uw,  vw,  ww  Wave  orbital  velocity  at  the  center  of  the  prop  disk,  defined 
in  the  ship-fixed  frame 

wp  Wake  fraction  coefficient 

xp,  yp,  zP  Position  of  the  center  of  the  propeller  disk  in  the  ship-fixed 
reference  frame. 

Z,(  Wave  elevation  above  propeller,  defined  in  earth-fixed  frame 

as  Angle  between  the  x-axis  of  the  ship-fixed  frame  and  the 

projection  of  the  shaft  axis  onto  the  vertical  plane  of 
symmetry  of  the  ship.  A  positive  angle  indicates  that  the 
shaft  points  downwards 

f]s  Angle  of  the  propeller  shaft  in  the  horizontal  plane  relative  to 

the  ship  x-axis,  positive  to  port. 

y  Flow  straightening  factor,  used  to  adjust  </>s  to  account  for  the 

flow  straightening  influence  of  the  hull. 

p  Density  of  the  fluid 

<j>s  Angle  between  shaft  axis  and  total  velocity  vector  at  the 

center  of  the  propeller  disk 
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2  Theory  Overview 

The  core  ‘physics  engine’  within  the  time-domain  simulator  solves  the  system  of  ordinary 
differential  equations  that  describe  the  equations  of  motion  (for  rigid  bodies)  and  propagates  the 
solutions  forward  in  time.  The  time-domain  simulator  and  equations  of  motion  are  described  in 
Section  3.  At  each  time  step  the  total  force  on  the  ship  must  be  computed.  In  order  to  compute 
the  forces  on  the  ship,  the  environment  in  which  the  ship  is  operating  must  first  be  evaluated, 
which  may  include  the  wind  and  current  in  addition  to  the  wave  environment.  The  environment 
model  will  be  able  to  compute  the  pressure  and  velocity  field  at  any  arbitrary  point.  In  the 
“Level-0”  implementation  of  Tempest,  the  environment  consists  only  of  the  superposition  of 
linear  long  crested  waves  as  described  in  Section  4.  When  the  ship  is  operating  in  waves  it  is 
necessary  to  control  the  rudders  to  keep  the  ship  on  course.  Tempest  uses  a  PID  controller  for 
the  rudder  which  is  described  in  Section  5.  In  Tempest  it  will  be  assumed  that  the  total  force  on 
the  ship  can  be  decomposed  into  various  force  contributions,  and  the  total  forces  on  the  hull  at 
each  time  step  can  therefore  be  obtained  by  summing  the  various  force  contributions.  In  the 
Level-0  theory,  the  following  force  contributions  will  be  considered: 

•  Hydrostatic  forces 

•  Froude-Krylov  forces 

•  Radiation  forces 

•  Diffraction  forces 

•  Hull  resistance 

•  Bilge  keel  forces 

•  Bare  hull  maneuvering  forces  (includes  skeg  force) 

•  Rudder  and  fin  forces 

•  Propulsion  forces 

•  Weight  of  ship 

The  total  force  on  the  ship  at  each  time  step  is  then  computed  as  the  sum  of  component 
forces.  The  theory  used  to  compute  each  component  force  is  described  in  Section  6.  The  total 
force  will  be  used  to  solve  the  rigid  body  equations  of  motion  in  ship-fixed  coordinate  system  at 
each  time  step. 
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3  Equations  of  Motion 


3. 1  Summary 

Tempest  will  use  a  time  stepping  method  to  solve  the  Euler  equations  of  motion  to  predict 
the  motion  of  a  ship.  This  section  presents  the  theory  for  setting  up  the  equations  of  motion. 
The  equations  of  motion  are  derived  from  Newton’s  Second  Law: 

F  =  md  (3.i) 

M  =  [I]a 

Equation  (3.1)  must  be  applied  in  an  inertial  coordinate  system.  It  would  be  possible  to 
solve  these  equations  directly  in  the  inertial  coordinate  system,  however  the  inertia  tensor,  [/],  as 
well  as  many  of  the  force  components,  is  defined  in  a  ship-fixed  axis  system  and  would  have  to 
be  derived  or  transfonned  into  the  inertial  system  at  each  time  step.  Therefore,  transformations 
are  applied  to  develop  the  equations  of  motion  in  a  ship-fixed  coordinate  system,  resulting  in  the 
Euler  equations  of  motion  for  a  rigid  body.  The  derivation  of  the  Euler  equation  of  motion  can 
be  found  in  several  sources  including  Abkowitz  (1969)  and  Etkin  (1982). 

The  force  and  moment  vectors  in  Equation  (3.1)  represent  the  total  force  and  moment 
acting  on  the  vessel.  In  this  document,  the  tenn  “forces”  will  hereafter  be  assumed  to  include 
both  forces  and  moments.  The  total  force  on  the  ship  at  each  time  step  will  be  computed  by 
summing  the  component  forces  described  in  Section  6. 

The  equations  will  first  be  formulated  making  no  assumptions  concerning  the  symmetry  of 
the  ship  or  the  location  of  the  center  of  gravity  of  the  ship  relative  to  the  origin.  In  the  Level-0 
theory,  however,  it  will  be  assumed  that  the  ship  is  symmetric  and  that  the  origin  of  the  ship- 
fixed  coordinate  system  coincides  with  the  center  of  gravity.  These  assumptions  will  be  used  to 
generate  a  simplified  set  of  equations  for  implementation  in  the  Level-0  Tempest  code.  No 
assumptions  are  made  that  the  ship-fixed  axis  system  is  aligned  with  the  principal  axes  of  the 
ship,  even  in  Level-0.  This  results  in  additional  tenns  containing  the  cross  products  of  inertia, 
which  are  often  assumed  to  be  zero. 

The  weight  of  the  ship  and  the  forces  from  the  infinite  frequency  added  mass  tenns  are 
included  in  the  equations  of  motion  presented  in  this  section.  The  theories  for  the  computation 
of  all  the  other  force  components  are  described  in  Section  6.  The  equations  described  in  this 
document  assume  that  the  forces  acting  on  the  ship  have  been  transformed  to  the  ship-fixed 
coordinate  system. 

Propagating  the  solution  to  the  equations  of  motion  is  the  job  of  a  numerical  integrator. 

This  propagation  is  often  done  using  a  linear  Euler  integration  scheme,  (i.e.  xM  =  xt  +  A tx  ), 

where  the  dot  over  the  x  indicates  differentiation  with  respect  to  time.  Though  this  simple 
scheme  is  adequate  for  non-oscillatory  systems  when  used  with  a  sufficiently  small  time  step,  in 
general,  higher  order  integrators  should  be  used.  Of  course,  more  accurate  integration  schemes 
tend  to  involve  more  calls  to  the  underlying  force  models  and  equations  of  motion.  Four 
different  integrators  were  developed  for  Tempest:  Euler,  Modified  Euler,  4th  Order  Runge-Kutta, 
and  Adams-Bashforth.  Additional  integrators  could  easily  be  added  as  new  components.  The 
equations  are  formulated  initially  for  the  Level-0  code,  with  the  intention  that  they  can  probably 
be  extended  with  only  minor  modifications  to  the  Level-II  and  Level-Ill  theory  implementations 
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of  Tempest.  Some  of  the  assumptions  regarding  symmetry  and  the  location  of  the  center  of 
gravity  may  be  removed  in  the  Level-II  and  Level-Ill  implementations. 


3.2  Coordinate  Systems  for  Equations  of  Motion 

The  equations  of  motion  are  set  up  to  describe  the  trajectory  of  the  ship-fixed  reference 
system  relative  to  an  earth-fixed  reference  system.  The  earth-fixed  reference  frame,  Ox  Y  z_ ,  is 

assumed  to  be  a  right-handed  coordinate  system  with  the  Z-axis  positive  upwards  and  the  origin 
on  the  calm  water  surface.  The  ship-fixed  reference  frame,  O  ,  is  a  right  handed  system  with 

the  x-axis  positive  forward  through  the  bow,  v-axis  positive  to  port,  and  the  z-axis  positive 
upwards.  In  the  Level-0  implementation  of  Tempest,  the  origin  of  the  ship-fixed  system 
coincides  with  the  center  of  gravity  of  the  ship.  Later  implementations  may  remove  this 
requirement,  so  the  equations  will  first  be  developed  for  an  arbitrary  location  of  the  origin.  The 
O  v  _  axis  system  is  fixed  with  the  ship  and  moves  with  all  the  motions  of  the  ship.  The 

0XeYe/r  axis  system  is  fixed  to  the  earth.  A  third  axis  system,  o  ,,  is  required  to  compute 

some  of  the  linearized  seakeeping  forces.  This  axis  system  moves  in  the  earth-fixed  X  and  Y 
direction  and  rotates  with  yaw  to  align  with  the  origin  of  the  ship-fixed  system.  This  third 
system’s  x-y  plane  remains  parallel  to  the  calm  water  surface,  and  it  does  not  heave,  pitch,  or  roll 
relative  to  the  earth-fixed  system.  Its  origin  is  the  same  as  the  ship-fixed  origin,  though  in  some 
special  cases  may  be  directly  below  the  ship-fixed  origin  on  the  calm  water  surface  (Ze  =  0). 

The  solution  to  the  equations  of  motion  will  provide  the  position  of  the  ship-fixed  reference 
frame  relative  to  the  earth-fixed  frame  at  each  time  step.  This  can  be  described  with  six 
variables:  the  ( X,Y,Z)  position  of  the  origin  of  the  ship-fixed  reference  frame  relative  to  the  earth- 
fixed  frame  and  the  three  Euler  angles:  (j),9  and  i//,  which  represent  the  roll,  pitch,  and  yaw  of  the 
ship-fixed  reference  frame.  The  order  of  the  rotations  are  important  and  must  be  performed  in 
the  order:  y/,  9,  (/).  When  the  rotations  are  taken  in  this  order,  the  definitions  of  yaw  ( y/),  pitch 
(9),  and  roll  ((f))  are  as  follows: 

•  yaw  ( y/)\  The  angle  between  the  vertical  plane  defined  by  the  earth- fixed  X,  Z  plane  and 
the  earth- fixed  vertical  plane  passing  through  the  ship-fixed  x  axis  (the  plane  defined  by 
the  ship-fixed  x  axis  and  the  earth-fixed  Z  axis  translated  to  the  ship-fixed  origin). 

Positive  yaw  is  defined  using  a  right  hand  rule  about  the  Z  axis. 

•  pitch  (9):  The  angle  of  elevation  of  the  ship-fixed  x  axis  relative  to  the  earth- fixed  X,  Y 
plane.  Positive  pitch  is  bow  down. 

•  roll  (</>):  The  angle  between  the  ship-fixed  x,z  plane  and  the  earth-fixed  vertical  plane 
passing  through  the  ship-fixed  x  axis  (the  plane  defined  by  the  ship-fixed  x  axis  and  the 
earth- fixed  Z  axis  translated  to  the  ship-fixed  origin).  Positive  roll  angle  is  defined  using 
a  right  hand  rule  about  the  x  axis.  Positive  roll  is  to  starboard. 

When  transforming  translational  velocities,  forces  and  moments  from  the  ship-fixed  frame 
to  the  earth- fixed  frame,  the  transformation  matrix  \Ts/e\  is  used,  where  this  matrix  is  defined  as: 


cos6*cos^ 
cos  9  sin  yr 
-sin  9 


-  cos^sin^  +  sin  9  sin  (f>  cos  yr  sin^sin^  +  sin  9  cos  (f)  cos  y/ 

cos  (f)  cos  y/  +  sin  9 sin  (f) sin  y/  -  sin^cos^  +  sin  9  cos (f>  sin  y/ 

cos  9  sin  (f)  cos  9  cos  (f) 


(3.2) 
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Alternatively,  the  quaternion  representation  can  be  used  to  fonnulate  the  position  and 
orientation  of  the  ship-fixed  frame  relative  to  the  earth-fixed  frame  can  be  formulated  in  terms  of 
the  (X,  Y,Z)  position  of  the  origin  of  the  ship-fixed  reference  frame  relative  to  the  earth- fixed 
frame  and  the  quaternion,  s.  The  quaternion  is  formed  by  four  components: 

s  \ 

6i 

s*  (3.3) 

£•3 

The  transfonnation  matrix  \Ts/e\  can  be  formulated  in  terms  of  the  quaternion  as: 

1  2(^£2  £3^  2(£,i£,2  £3£4)  2(^£ x£ 3  ~\~  £ 2£ 4^)  (3  4) 

S  /  £■  ]  ^(^1^2  ^3^4  )  i  —  2(f|  +£3)  2{£2£3  —  ) 

2(£l£3-£2£4)  2(£2£3+£l£4)  l~2(£2+£2) 


The  quaternion  components  satisfy  the  constraint  equation: 

Sx  +  £  2  +  £3  +  £  4  =  1  (3 

The  quaternion  components  can  be  computed  from  the  Euler  angles  and  vice  versa.  The 
equations  to  obtain  the  quaternion  from  the  Euler  angles  are: 


For  most  cases  the  initial  roll,  pitch  and  yaw  will  be  zero,  in  which  case  the  initial  value  of  the 
quaternion  is  (0,0,0, 1).  The  equations  used  to  obtain  the  Euler  angles  from  the  quaternion  are: 


(j)  =  arctan 


2(g2g3  +g,g4)" 

l-liel+el), 


0  =  arcsin (2(£2£4  -£l£i)) 


y/  =  arctan 


2(g,g2  +g3g4)^ 

l-2(s22  +  s32). 


(3.7) 
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The  transformation  matrix  [Ts/e]  is  used  to  transform  the  position  or  translational  velocity  of 
point  in  the  ship-fixed  frame  to  the  earth- fixed  frame.  To  compute  the  earth- fixed  translational 
velocities  the  formula  is: 


rx\ 

(u \ 

Y 

-  ITs/^] 

V 

(3.8) 


To  convert  the  x,y,z  position  of  a  point  defined  in  the  ship-fixed  frame  to  the  earth-fixed  frame, 
the  formula  is: 


Y 

Z 

v  /Earth -fixed 


Y 


+  1 Ts/e  ] 


y 


v  J  Ship -fixed 


(3.9) 


where  the  (. X,  Y,Z)  vector  on  the  left  hand  side  is  the  position  of  the  point  in  the  earth-fixed  frame, 
the  first  ( X \  Y,Z)  vector  on  the  right  hand  side  is  the  position  of  the  origin  of  the  ship-fixed  frame 
defined  in  the  earth-fixed  frame,  and  the  last  ( x,y,z )  vector  is  the  position  of  the  point  in  the  ship- 
fixed  frame.  The  matrix  [  Ts/e]  is  orthogonal,  so  its  inverse  is  simply  the  transpose  of  the  matrix, 
which  can  be  used  to  perform  transformations  in  the  opposite  direction  to  convert  quantities  from 
the  earth-fixed  frame  to  the  ship-fixed  frame. 


cos  #  cos  y/  cos  #  sin  y/  -sin# 

-  cos  (j)  sin  y/  +  sin  #  sin  (j)  cos  y/  cos  (j)  cos  y/  +  sin  #  sin  (j)  sin  y/  cos#  sin  </> 
sin^sin^  +  sin#  cos  ^  cos  ^  -  sin  (j)  cos  y/  +  sin  #  cos  (j)  sin  y/  cos#cos  (j) 


1  2(£’2  +£3)  2(£-j£-2  +  £'3 £  A  ) 

2 ( /.’ | E 2  6*2^4)  1  2(6*^  6*2  ) 

2(^1  ^3  +^2^4)  2(^2^3-^1^4) 


2(^1  ^3  -^2^4) 
2(£2£3  +£l£4) 
1  -  2(&'f  +^2) 


and: 


M 

V 

=fr„] 

Y 

A 

(3.10) 


(3.11) 


A  different  set  of  equations  are  required  to  compute  the  time  derivatives  of  the  Euler  angles 
from  the  angular  velocities  defined  in  the  ship-fixed  axis  system  and  vise-versa.  Note  that  the 
derivatives  of  the  Euler  angles  are  not  the  rotational  velocities  with  respect  to  the  earth-fixed 
reference  frame,  but  rather  the  rate  of  change  in  the  roll,  pitch  and  yaw  as  defined  above.  The 
Euler  angle  derivatives  do  not  form  an  angular  velocity  vector,  as  they  do  not  represent  rotations 
about  an  orthogonal  set  of  axes.  Therefore  the  matrix  that  is  used  to  transform  is  not  orthogonal. 
The  following  transfonnation  matrix  is  used  to  compute  the  rotational  velocity  vector  in  the  ship- 
fixed  frame,  (p,q,r ),  when  the  Euler  angle  derivatives  (<j>,  #,  ///) ,  are  known: 
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(3.12) 


V 

“1 

0 

q 

= 

0 

cos^ 

0 

-  sin^i 

-  sin  9 
cos  6  sin  (f> 
cos  6  cos  (f> 


f  )\ 

<t> 

0 


The  matrix  can  be  inverted  to  compute  the  Euler  angle  derivatives  (<j>,  0,y/  )  from  the 
rotational  velocity  vector  in  the  ship-fixed  frame  ( p,q,r ): 


(*} 

1 

tan  9  sin  (j) 

tan  9  cos  (j) 

9 

= 

0 

COS  (/) 

-  sin  (j) 

q 

j) 

0 

sin  </>/ 

/ cos  9 

COS  (j)/ 

/ cos  9  _ 

There  is  a  singularity  in  the  transformation  matrix  shown  above  when  the  pitch  angle  of  the 
ship  is  ±90°,  however  this  is  not  expected  to  be  a  problem  for  a  surface  ship,  even  if  the  ship  is 
operating  in  very  large  waves.  This  problem  can  be  avoided  when  the  quaternion  fonnulation  is 
used.  With  the  quaternion  representation,  the  time  derivative  of  the  quaternion  is  computed 
instead  of  the  time  derivative  of  the  Euler  angles.  The  time  derivate  of  the  Euler  angles  can  be 
computed  from  ( p,q,r ): 


(3.14) 


Finally  a  transformation  matrix  may  be  required  to  transform  forces  and  velocities  from  the 
yawed  upright  axis  system,  0rV_,,  to  the  ship-fixed  axis  system.  The  matrix  used  to  transform 

quantities  from  the  yawed  upright  frame  to  the  ship-fixed  frame  is  the  matrix  \Tu/s\  which  is 
shown  below.  This  matrix  is  orthogonal,  so  its  inverse  is  simply  the  transpose  of  the  matrix. 


cos  9 
sin  <f>  sin  9 
cos  (j)  sin  9 


0  -  sin  9 

cos  </>  sin  (f)  cos  9 
-sin^  cos  (j)  cos  9 


(3.15) 


[Ts,v]  =  {Tuls]T 


cos  <9 

sin  (j)  sin  9 

cos  (f)  sin  9 

0 

COS  (f) 

-sin^ 

-sin  9 

sin  (j)  cos  9 

cos^cos  9 

(3.16) 
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Figure  3-1  Definition  of  reference  frames  used  in  equations  of  motion 


3.2.1  Heading  Definitions 

Tempest  allows  the  ship  to  head  in  any  direction  within  the  earth-fixed  coordinate  system, 
and  likewise  allows  waves  to  travel  in  any  direction.  The  ship  heading  is  defined  as  the  angle  y/, 
which  is  equivalent  to  the  yaw  angle  defined  in  the  Section  3.2.  The  heading  of  a  wave  system 
in  earth-fixed  coordinates  is  defined  by  the  angle  fl  The  relative  wave  heading,  which  is  the 
heading  of  the  wave  system  relative  to  the  ship’s  heading,  is  given  by: 

Mn  =  Pn~¥  (3-17) 

The  ship  heading,  wave  heading,  and  relative  wave  heading,  are  shown  in  Figure  3-2. 
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Note:  Z  is  +  out  of  the  page 


y/=  Ship  heading 
Po  -  Primary  Wave  heading 
pn  =  Wave  Heading  of  nth 
wave  system 

Ho  =  Relative  wave  heading 

V - - S 

“Heading  of  primary 
wave  relative  to  ship” 

Not  Shown: 

jun  =  Relative  nth  wave  heading 

fin  =  Pn  -  v 

J3n,o  =  nth  wave  system  heading 
relative  to  primary  waves 

fti.O  —  Pn~  Po 


Figure  3-2  Definition  of  ship  heading  and  wave  heading  angles 


3.3  General  Formulation 

In  Tempest  the  equations  of  motion  using  both  an  Euler  angle  formulation  and  a  quaternion 
formulation  will  be  implemented.  The  majority  of  the  implantation  is  identical  for  both 
formulations,  with  the  primary  difference  being  that  the  time  integration  is  performed  on  the 
derivatives  of  the  Euler  angles  or  quaternion  components.  The  equations  of  motion  for  a  six- 
degree-of- freedom  rigid  body  can  be  derived  from  Newton's  Laws  of  Motion  (see  Abkowitz 
(1969)  and  Etkin  (1982),  for  example).  If  the  mass  of  the  body  and  the  mass  distribution  are 
constant  over  time,  then  the  equations  can  be  written  as: 

m  (u  +  qw-  rv -  xG(q 2  +  r2)  +  yG(pq -r)+ zG(pr +q))=Fx  +  W sin# 


m  ( v+ru-pw-yG(r 2  +  p2)  +  zG (qr - p)  +  xG(qp+ r)]=  Fy  -fEcos#sin^ 


m(w+ pv-qu-zG(p2  +q2)  +  xG(rp-q)  +  yG(rq  +  p])  =  F_  -  WcosGcosrj) 

I xx  P  +  if  zz  -  i yy) qr  +  (  P>'  ~  d)  I xy  ~  {f+pq)lzx  +  (r2  ~  q^)  I yz  +  « [»,  (w  ~  W  +  vp)  ~zG(v~wp+  ur )] 

=  Mx  -  y CW  cos  6  cos  (f>  +  zcW  cos  6  sin  (j) 

lyy  q  +  {ixx  -  I  zz)rp  -  (  p  +qr)  I  w  +  ip2  -  r2)  I  zx  +  {qp-r)  I  r  +  m[zG  (it  -  vr  +  wq )  -xG(w-  uq  +  vp)] 
=  M r  +  xGW  cos  9  cos  (j)  +  zGW  sin  0 


17 


Izz  r  +  (fyy  -  Ixx)  PA  +  (  A  ~  I xy  +  {>P  ~  P  )  I zx  ~  (p  +  ’P)  I yz  +  m[XG  (v  ~  Wp  +  Ur)  -  (it  ~  VT  +  Wq )] 
=  M z  -  xcW  cos  #  sin  (f>  -  yGW  sin  # 


(3.18) 

The  Euler  angle  terms  that  appear  in  the  above  equations  result  from  the  transformation  of 
the  weight  force  vector  from  the  earth-fixed  frame  to  the  ship-fixed  frame.  If  the  quaternion 
formulation  is  used,  the  relations  shown  below  can  be  substituted  to  fonnulate  the  equations  in 
terms  of  the  quaternion  components. 


sin#  =  2 (s2s4  - £, s3 ) 
cos#  sin  (j>  =  2(s2£i  +  £x£4) 
cos  #  cos  (j)  =  1  -  2(p2  +  £p) 


(3.19) 


It  has  been  found  that  moving  the  terms  related  to  the  infinite  frequency  added  mass  to  the 
left  hand  side  of  the  equations  of  motion  improves  the  stability  of  the  system.  This  will  be 
implemented  in  Tempest.  In  the  Level-0  implementation  it  will  be  the  total  added  mass, 
obtained  by  interpolating  the  added  mass  at  the  wave  encounter  frequency,  that  will  be  included 
on  the  left  hand  side  of  the  equations  of  motion,  instead  of  the  infinite  frequency  added  mass 
terms.  In  almost  all  cases,  the  geometry  of  the  ship  will  be  symmetrical,  in  which  case  many  of 
the  added  mass  terms  in  the  above  matrix  will  be  zero.  Any  added  mass  coefficient  with  one 
even  subscript  and  one  odd  subscript  will  be  zero  for  an  upright  symmetric  ship  on  the  calm 
water  surface.  For  the  “Level-0”  theory,  the  added  mass  and  damping  coefficients  will  be 
computed  using  a  method  that  is  linearized  about  this  condition,  and  in  this  case  A /_?.  A I4,  A24,  etc 
can  be  assumed  to  be  zero.  However,  in  higher  levels  of  theory  this  linear  approximation  may 
not  be  used,  and  symmetry  cannot  be  used  when  computing  the  added  mass  and  damping  terms 
for  a  heeled  ship.  For  this  reason,  all  the  terms  in  the  added  mass  and  damping  matrices  will  be 
included  in  the  general  derivation.  A  simplified  derivation  for  implementation  in  “Level-0”  will 
be  provided  in  the  following  section. 

The  six  equations  will  be  formulated  to  solve  for  the  acceleration  terms  in  each  mode,  with 
terms  related  to  the  velocity  and  displacement  in  each  mode  appearing  on  the  right  hand  side, 
computed  using  values  from  previous  or  intermediate  time  steps  depending  on  the  time 
integration  scheme  used.  In  matrix  form  the  equations  will  appear  as  shown  below: 


a\l 

*12 

*13 

*14 

*15 

*16 

U 

1 

"2 

l _ 

a2\ 

*22 

*23 

*24 

*25 

*26 

V 

rhs2 

*31 

*32 

*33 

*34 

*35 

*36 

w 

rhs. 

a4\ 

^42 

*43 

*44 

*45 

*46 

p 

rhs4 

*51 

*52 

*53 

*54 

*55 

*56 

q 

i'hs  5 

_«61 

*62 

*63 

*64 

*65 

*66  _ 

r 

rhs6 

(3.20) 
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The  coefficients  of  the  left  hand  side  matrix  are  defined  below. 


an  =  m+Auoo 
an  —  A]2oo 
013  =  Ai3oo 
a  14  =  A]4oo 

a  15  =  zGm  +A1500 
a  i6  =  -  yam  +AJ6co 

aii  =  A2100 
a2i  =  m  +  A2200 
a  23  =  A2300 

a24  ~—  ( i  01  A2400 

025  =  A2500 

Cl26  ~  Xfjm  +  A2600 

031  =  A2100 
032  =  A3200 
a33  =  m  +  Ajsoo 
d34  =  VG  m  +  A3400 

CI35  =  -Xg  m  +A3500 

a36  A3600 

041  A4I00 

a  42  =  -zGm  +  A 42  co 
043  =  ycm  +  A43oo 

a44  f 44  A4400 

O45  ~Ixy  A4500 

O46  Izx  A4600 

051  =  ZGm  +  As  ioo 
052  =  A5200 
053  =  -xGm  +  A5300 
O54  Ixy  T  A5400 

055  lyv  A5500 

O56  Izx  T  A5600 

061  =  -vg  01  +  A  61  oo 

062  -  XGm  +  A62oo 

063  A  63  oo 

064  ~Izx  A64oo 

065  ~Iyz  A65oo 

066  Izz  T  A  6600 
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rhsx 

rhs2 

rhs3 

rhs4 

rhs5 

rhs6 


=  F'x  +  W  sm0-mqw-rv-%(q2  +r2  )+yGpq+zGpr\ 

=  F'v  -  W cosO s'm (f>- m  ru- pw-yG(r2  +  p2 )  +  zGqr  +  xGqp\ 

=  F'z  -  WcosGsmG-m  pv-qu-zG(p1  +q2  )+xGrp+yGrq) 

=  M'x-ycw  cos  9 COS (/)  +  zGW  cos  9 sin  </>  -  (/_,  -  /,,,,) qr  -  pr  Jxy  -  pq  Izx  -  (r2  -  q2)  Iyz 

-  m  \yG  (; uq  +  vp )  -  zG  ( wp  +  ur )] 

=  M'y  +  xGW  cos  9  cos  (p  +  zGW  sin  0  -  (fxx  -  Jz:)  rp  -  qr  fxy  -  [p2  -  r2)jzx  -  qp  Jy: 
-m\zG{vr+  wq) -xG(uq  +  vp)\ 

^M'z-xGWcos9sm</)-yGWsm9-  {lyy  -  fxx) pq  -  ( q2  -  p2)lxy  -  rq  fzx  +  (rp)  [yz 

-  m[xG  (m p  +  ur)  -  yG  (vr  +  wq)\ 


(3.21) 


F'  and  M'  are  the  total  forces  and  moments  not  including  the  weight  of  the  ship  and  the 
force  contribution  from  the  infinite  frequency  added  mass  terms. 

At  every  time  step,  the  six  equations  listed  above  in  matrix  form  are  solved  to  compute  the 
accelerations  in  the  ship-fixed  frame,  (u,  v,  w,  p,  q,  r)  .  The  ship-fixed  accelerations  are  integrated 
with  time  to  obtain  the  ship-fixed  velocities  ( u,v,w,p,q,r ). 

u  =  u0+^udt  ,  v  =  v0  +  Jvr/f  ,  w  =  wQ+^wdt  (3.22) 

P  =  Po+\pdt  ,  q  =  q0+jqdt  ,  r  =  r0+^rdt 


Here  (: uo,vo,wo,Po,ro,qo )  are  the  initial  translational  and  rotational  velocities  in  the  ship-fixed 
frame. 

Up  to  this  point  the  procedure  is  identical  for  the  Euler  angle  formulation  and  the 
quaternion  formulation.  The  next  step  is  where  the  two  formulations  differ.  In  the  Euler  angle 
formulation,  the  time  derivatives  of  the  Euler  angles  derivatives  (<f>,9,y/)  are  computed  from  the 
ship-fixed  rotation  velocities  (p,q,r ),  using  Eqn.  (3.13).  The  time  integration  is  perfonned  on  the 
Euler  angle  derivatives  to  obtain  the  Euler  angles  at  the  current  time  step: 

0  =  (f>o+^(f>dt  ,  0  =  Gq  +^Gdt  ,  y/  =  y/0+^ipdt  (3.23) 


where  <f>o,  Go,  and  y/o  are  the  initial  roll,  pitch  and  yaw  angles  of  the  ship. 

In  the  quaternion  formulation,  the  time  derivatives  of  the  quaternion  components  are 
computed  instead  of  the  time  derivatives  of  the  Euler  angles,  and  then  the  time  integration  is 
performed  for  the  quaternion  instead  of  the  Euler  angles. 
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(3.24) 


£\  —  £y  0  +  j*  Gy  dt  ,  G  T  —  £2  Q  +  J  £~)  dt 

£3  =  £3  0  4"  J"  £3  dt  ,  £4  =  £4  0  +  j  £4  dt 

The  _0  subscript  indicates  the  initial  value.  For  the  quaternion  formulation,  numerical 
integration  errors  accumulate  so  that  the  constraint  imposed  by  in  Eqn.  (3.5)  is  no  longer 
satisfied.  The  quaternion  should  be  re-normalized  at  regular  intervals  during  the  simulation  to 
satisfy  Eqn.  (3.5). 

The  transformation  matrix  \Ts/e\  is  then  computed  using  either  Eqn.  (3.2)  in  the  Euler  angle 
formulation  or  Eqn.  (3.4)  in  the  quaternion  fonnulation.  The  ship-fixed  translational  velocities 
( u,v,w )  are  then  transformed  to  the  translational  velocities  in  the  earth-fixed  frame  ( Sc,y,z ) ,  and 
then  integrated  to  obtain  the  position  of  the  ship  relative  to  the  earth-fixed  frame. 

X  =  X0+jxdt  ,  Y  =  Y0  +  jfdt  ,  Z  =  Z0+jzdt  (3  25) 

where  (Xo,Yo,Zo)  is  the  initial  position  of  the  ship-fixed  origin  defined  in  the  earth- fixed  frame. 

3.4  Restricting  degrees  of  freedom 

In  certain  cases  it  may  be  desired  to  restrict  the  ship  in  one  or  more  degrees  of  freedom. 
This  may  be  necessary  to  simulate  a  model  test  where  the  model  was  free  to  move  only  in  some 
degrees  of  freedom,  or  to  compute  forces  on  a  ship  with  a  prescribed  motion.  The  restriction  on 
the  motion  of  the  ship  can  be  applied  either  in  the  earth-fixed  frame  or  in  the  ship-fixed  frame. 
In  the  case  where  the  restrictions  are  imposed  in  the  ship-fixed  frame,  the  time  history  of  the 
ship-fixed  acceleration  in  the  restricted  degrees  of  freedom  will  be  provided  as  input. 
Alternatively  the  time  history  of  the  ship-fixed  velocity  can  be  specified  as  input  and  the 
accelerations  computed  by  differentiation.  Eqn.  (3.20)  will  still  be  used  to  solve  for  all  six 
accelerations  in  the  ship-fixed  frame.  However,  the  ship-fixed  acceleration  component 
corresponding  to  restricted  degrees  of  freedoms  u,v,w,p,q  or  r  will  be  set  to  the  user-specified 
value  prior  to  performing  the  integration  to  compute  the  ship-fixed  velocities.  The  computed 
values  for  ship-fixed  accelerations  from  solving  Eqn.  (3.20)  can  be  used  to  compute  the  forces 
and  moments  that  must  be  applied  to  restrict  the  ship  in  the  specified  degree  of  freedoms. 

In  some  instances,  such  as  simulating  a  model  test  where  the  ship  is  restricted  by  using  a 
heave  staff,  it  is  the  motion  in  the  earth-fixed  frame  that  is  restricted.  In  this  case  one  or  more 
components  of  translational  velocity  vector  in  the  earth-fixed  frame,  ( X,Y,Z ) ,  or  one  or  more  of 
the  time  derivatives  of  the  Euler  angles,  ,  must  be  provided  as  input.  For  this  case,  the 

ship-fixed  accelerations  and  velocities  will  first  be  computed  without  any  restrictions.  Then  the 
unrestricted  values  for  (X,Y,Z)  are  computed  using  Eqn.  (3.8)  and  the  unrestricted  values  for 
(<j>,0,y/)  are  computed  using  Eqn.  (3.13).  The  values  for  (X, Y, Z)  and  (cj),9,y/)  in  the  restricted 
modes  are  then  replaced  with  their  prescribed  values.  The  ship-fixed  velocities,  ( u,v,w,p,q,r ), 
including  the  effects  of  the  restricted  motion  can  then  be  computed  using  Eqns.  (3.11)  and  (3.12). 

3.5  Simplified  Formulation  for  Level-0  Implementation 

In  the  Level-0  implementation  of  Tempest,  several  assumptions  will  be  made  which  result 
in  a  simpler  form  of  the  equation  of  motion.  These  assumptions  are: 


21 


The  origin  of  the  ship-fixed  frame  is  at  the  center  of  gravity  of  the  ship 

The  geometry  of  the  ship  is  symmetric,  so  any  added  mass  tenns  with  both  an  even  and 
odd  subscript  can  be  set  to  zero. 

The  complete  added  mass  interpolated  at  the  encounter  frequency  will  be  included  on  the 
left  hand  side  of  the  equations  in  Level-0,  as  described  in  the  Section  6.2. 1.3. 


With  these  assumptions,  the  equations  of  motion  simplify  to: 

m  (u  +  qw-rv)=  Fx  +  W  sin  6 
m  (v  +  ru  -  pw )  =  Fv  -W  cos  9  sin  <j> 
m  (w  +  pv-  qu)  =  Fz-W  cos  6  cos  <f) 

I xx  P  +  (izz  -  Iyy)  <7^  +  (pr~q)  I xy  ~  0"+m)  I zx  +  (/  ~  S)  I  yz  =  ^  x 
Iyy  q  +  (ixx  -  Izz)  rP~{p  +  P1)  I  xy  +  (/?'  "  I  zx  +  (pP  ~  T  )  I  yz  =  Mv 

Izz  r  +  (iyy  -  Ixx) pq  +  (q2-  p  )  I xy  +  (rq  -  p  )  I zx  -  (q  +  rp)  J yz  =  Mz  (3.26) 


Adding  the  added  mass  terms  described  in  Section  6. 2. 1.3  to  the  left  hand  side  and 
rewriting  in  matrix  form,  these  equations  can  be  written  as: 
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where, 


(3.27) 


rhs ,  =  Fx  '+W  sin  6  -  m[qw  -  rv] 

rhs2  =  Fy  '-W  cos 6 sin (/)  -  m[ru  -  pw] 

rhs  3  =  F_  '-W  cos  9  sin  9  -  m[pv  -  qu  ] 

rhs,  =  Mx '-(/zz  - 1^ )qr  -  prlxv  -  pqlzx  -  (r2  -  q2  )Iyz ) 

rhs 5  =  M;-(4  -  Izz)rp  -  qrlxy  -  (p1  -  r2)Izx  -  qplyz) 

rhs6  =  M: '-(i^  -  4 )pq  —  (q2  —  p2  )lxy  -  rqlzx  -  rplyz 
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F'  and  M'  are  now  the  total  forces  and  moments  not  including  the  weight  of  the  ship  and  the 
force  contribution  from  the  added  mass  terms.  F'  and  M'  are  computed  by  summing  the 
component  forces  described  in  Section  6  according  to  Equation  (6.2).  The  total  added  mass 
force  contribution  is  now  included  on  the  left  hand  side  by  the  Ay  ’  terms.  After  computing  the 
accelerations  in  the  ship-fixed  frame  (u,v,w,p,q,r) ,  the  procedure  for  integrating  the 
accelerations  to  compute  the  velocities  and  integrating  the  velocities  to  compute  the  ship 
trajectory  is  identical  to  that  listed  in  the  previous  section. 

4  Linear  Wave  Environment  for  Tempest 


4.1  Introduction 

This  section  details  the  theory  for  computing  the  wave  elevation,  pressures,  and  kinematics 
for  a  linear  seaway  to  be  used  in  the  Tempest  dynamic  stability  and  seakeeping  simulation  code. 
The  kinematics  for  a  linear  regular  (sine)  wave  travelling  in  an  arbitrary  direction  in  infinite  and 
finite  depth  water  are  presented.  The  definition  for  an  irregular  sea  is  built  upon  these 
kinematics  and  includes  both  long-crested  (unidirectional)  and  short-crested  (directional  spread) 
seaways.  The  formulations  for  several  theoretical  spectra  are  included.  Specifications  for  input 
and  output  parameters  follow  the  development  of  the  theory. 

4.2  Coordinate  System 

The  coordinate  system  used  to  describe  the  wave  field  is  an  earth-fixed  coordinate  system 
defined  by  a  right  hand  rule  with  Z  positive  up.  Z  is  equal  to  zero  on  the  calm  water  line.  A 
wave  direction  of  /?=  0  means  that  a  wave  progresses  along  the  X-axis. 


4.3  Relationships  for  Infinite  Depth  Linear  Regular  Waves 

Wave  Number: 


Wave  Length: 


Frequency: 


k 


2  -JT  _  CO2 
2  8 


2  •  7T  _  2  •  TZ  •  g  _  g  '  T~ 

k  co2  2  •  tc 


O  f  I — T  2'n  \2'n-g 

co  =  2-n-f  =  ^g-k=  —  =^—j  — 


Wave  Celerity  (Phase  Velocity): 

_  oo 
k 


g 

2-n-f 


Group  Velocity: 


(4.1) 

(4.2) 

(4.3) 

(4.4) 

(4.5) 
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4.4  Relationships  for  Finite  Depth  Linear  Regular  Waves 

Wave  Number: 


k  =  — —  =  k  ■  tanh(£  •  h) 

S 


Wave  Length: 


2-n  _2-n  ■  g  .  ,(  2  -n-h^ 


■  tanh 


CO 


g  •  T2  ■  tanh 


r  2  ■  n  ■  h  ^ 


\  A  J 


2- 


K 


\  A  J 


Frequency: 


co  =  2-n  ■  f  =  —  —  =  ^/g  •  k  ■  tanh(A:  •  h)  = 


Wave  Celerity  (Phase  Velocity): 


2-x-g 

X 


tanh 


2-n  ■  h 

V  A  J) 


c  =  —  =  J—  ■  tanh(£  •  h)  =  —  ■  tanh(A:  •  h )  =  — — — 
k  V  k  co  2-n-  f 


Group  Velocity: 


_  co  f  1  k  -h  ' 
8  k  v2  sinh(2  -k  -h)y 
c g  — >  c  as  A  — »  0 

1 

— >  —  c  as  /z  — »  oo 
g  2 


(4.6) 


(4.7) 


(4.8) 


(4.9) 


(4.10) 


4.5  Regular-Linear-Wave  Velocity  Potential  and  Kinematics 

The  3D  velocity  potential  of  a  regular  wave  traveling  in  an  arbitrary  direction  in  finite  and 
infinite  depth  is  presented  in  the  following  sections.  One  can  easily  reduce  these  equations  to  a 
2D  wave  by  setting  the  wave  direction,  /?,  to  0. 


4.6  3D  Velocity  Potential  and  Kinematics  for  a  Single  Wave  -  Infinite- 
Depth 

Velocity  Potential: 


q>  =  LAXL .  gk-z 

CO 


gi-(co-t-k-X -cos{p)-k-Y  -sin(/?)+£) 


(4.11) 


Real  Part  of  the  Velocity  Potential: 


24 


n<p)= 


(4.12) 


^  •  ekz  ■  sin(<z>  -t-k-X  ■  cos  (/?)  -k-Y  ■  sin(/?)  +  s ) 

co 


Wave  Elevation: 


c 


-1 

g 


( d  } 
x*(<P) 
\dt 


a  ■  cos(a>  -t-k-X  ■  cos(/?)  -k-Y  ■  sin(/?)  +  s) 


Hydrostatic  Pressure: 


Dynamic  Pressure: 

( d  ^ 

Pd=~P ■ 

\dt 


Phys=-p-g-Z 

p  ■  g  -  a-  ekz  ■  cos(a»  -t-k-X  ■  cos(/?)  -  k-Y  ■  sin(/?)  +  s) 


Particle  Orbital  Velocity  Components: 

F  =  {«w,vw,ww}=VSR(^)  =  f^ nitpxJLxitpxl-Viitp) 

\OJC  uY  uZ  y 

uw  =  a-  co-  ekZ  ■  cos(u)  -t-k-X  ■  cos(/?)  -k-Y  ■  sin(/?)  +  s )  •  cos(/?) 
vw  =  a-  co-  ekZ  ■  cos  (co  -t-k-X  ■  cos(/?)  -  k-  Y  ■  sin(/?)  +  s )  •  sin(/7) 
ww  =  -a-  co-  ekZ  ■  sin(n>  -t-k-X  ■  cos(/?)  -k-Y  ■  sin(/?)  +  s ) 
Particle  Acceleration  Components: 


ot  ot 

uw  =  -a  ■  co2  -ekZ  ■  sin(<y  •  t  -  k  ■  X  ■  cos  (/?)  -k-Y  ■  sin(/?)  +  f)  •  cos(/?) 
vw  =  a  ■  co1  ■  ekz  ■  sin(<y  •  t  -  k  ■  X  ■  cos  (J3)  -k-Y  -  sin(/?)  +  s)  -  si  n(j3) 
ww  =  -a  ■  co2  ■  ekZ  ■  cos(<y  -t-k-X  ■  cos  (/?)  -  k-  Y  ■  sin(/?)  +  s ) 


(4.13) 


(4.14) 

(4.15) 


(4.16) 


(4.17) 


4.7  3D  Velocity  Potential  and  Kinematics  for  a  Single  Wave  -  Finite 
Depth 

Velocity  Potential: 


i  ‘  g  '  ®  C0sh(  k  '  ( h  ~\~  Zy)  ^i-(co-t-k-X -cos(/?)- k-Y -sin(/?)+£) 

co  cosh(£  •  h) 


(4.18) 


Real  Part  of  the  Velocity  Potential: 

iR(^)  =  — — —  ■  C0S^^  ^  +  Z^  .  sin((y  -t-k-X  ■  cos  (/?)  -k-Y  ■  sin(/?)  +  s) 
co  cosh(A:  ■  h) 


(4.19) 


Wave  Elevation: 
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-1  ( d  \ 

C,  = - —  iR(( 'p )  =  a  ■  cos(&>  -t-k-  X  ■  cos (/?)  -k  -Y  ■  sin (/?)  +  s )  (4.20) 

g  vdt 

Hydrostatic  Pressure: 

P hys  =  P  '  g  '  Z  (4.21) 

Dynamic  Pressure: 

pD  = -p  ■( =  p-g-a-  cos^^  — + — ^  ■  cos(a>  -t-k  ■  X  ■  cos  (/?)  -k-Y  ■  sin  (/?)  +  g) 

D  U*  J  coshfiL//)  (4.22) 


Particle  Orbital  Velocity  Components: 

V  =  {uw,vw,ww}=  V5R(p)  =  f^*(P).^*(P).^tt(P) 
..  _a'£^  COSh(M/7  +  Z))  _  .  ,  v  ,,  v 


co  cosh  (k-h) 

a-g-k  cosh (k-Ui  +  ZY) 

v  =  — - - - — - - -  •  cos( 

co  cosh  (k-h) 

-a-g-k  sinh (k-Ui  +  ZY)  . 

Ww  = - 5 - —  •  SI 

co  cosh(k  ■  h ) 

Particle  Acceleration  Components: 


cos  {co  -t-k-  X  ■  cos(/?)  -k-Y  ■  sin(/?)  +  s )  ■  cos(/?) 


•  cos(<u  -t-k  ■  X  ■  cos  (/?)  -k-Y  ■  sin(/?)  +  s)  ■  sin(/?) 


•  sin(<u  -t-k  ■  X  ■  cos  (/?)  -k-Y  ■  sin(/?)  +  s) 


(4.23) 


V  =  {»>».  w„.}  =  — =  — (v«(^)) 

ot  at 

uw=-a-  g-k-  cosh(^  (h  +  Z))  s-n^ .  f  _  £ .  x  ■  cos(/?)  -k-Y  ■  sin  (/?)  +  £)  •  cos(/?) 
cosh(k-/7) 

v  =-a-  g-k  -  C0S^^  —  +  ^  •  sin(cu  -t-k-X  ■  cos(/?)  -k-Y  ■  sin  (/?)  +  s)  ■  sin(/?) 
cosh (k  ■  h) 

ww=-a-  g-k-  —  +  ^  ■  cos(co  -t-k-X  ■  cos  (/?)  -k-Y  ■  sin  (/?)  +  <?) 

cosh(k-/7) 

Note  that  the  finite-depth  dispersion  relation  is  different  than  the  infinite-depth  case  and 
one  cannot  substitute  k  =  co  /g  when  deriving/defining  the  kinematics  for  the  finite-depth  case. 


(4.24) 


4.8  Pressure  Stretching 

Linear  wave  theory  has  a  deficiency  in  that  the  total  pressure  (dynamic  plus  hydrostatic) 
does  not  approach  zero  as  the  free  surface  is  approached  from  below.  There  are  several  methods 
to  deal  with  this  issue  that  modify  the  dynamic  pressure,  orbital  velocity,  and  acceleration 
components. 

Wheeler's  Method  for  stretching  can  be  written  as  a  modification  to  Z  in  the  dynamic 
pressure  and  orbital  velocity  equations.  The  value  Z'  is  computed  and  used  in  place  of  Z  to 
compute  the  dynamic  pressure,  particle  velocities,  and  particle  accelerations.  The  hydrostatic 
pressure  is  not  stretched  and  uses  Z,  not  Z'. 
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for  infinite  depth 


(4.25) 


Z'=a(Z-£) 


a 


1 

h 


{ h  +  £ 


for  finite  depth 


In  the  case  of  irregular  waves,  C,  in  the  above  equation  is  the  total  free  surface  elevation 
comprised  of  all  wave  components  from  all  seaways  present  in  the  simulation. 


4.9  Irregular  Wave  Kinematics  -  Linear  Model 

Irregular  waves  are  computed  by  assuming  a  power  spectrum  for  a  seaway.  From  the 
spectrum,  individual  wave  components  are  generated.  The  characteristics  of  these  wave 
components  are  summed  to  compute  the  total  wave  characteristics. 

4.1 0 Long-Crested  Infinite  Depth  Irregular  Seaway  Kinematics 

An  irregular  seaway  is  modeled  as  a  superposition  of  regular,  sinusoidal  waves,  as  defined 
in  previous  sections.  The  amplitudes  of  the  wave  components  that  compose  the  irregular  seaway 
are  derived  from  a  power  spectrum,  S(co),  which  can  be  derived  from  measurements  or 
theoretical  spectra  (which  are  defined  in  a  subsequent  section). 

The  amplitudes  for  the  wave  components  are  derived  from  a  discretized  power  spectrum 
and  are  computed  as  follows: 

ai  =  -yj 2  ■  S(coi )  •  8oj! 

The  wave  elevation  at  an  arbitrary  location  is  the  computed  as 

Ctom  = 

i= 1 

The  phase  angles  (<q)  for  each  of  the  wave  components  are  random  numbers  which  are 
uniformly  distributed  between  (-n,n)  and  drawn  independently  for  each  frequency 

The  hydrostatic  pressure  is  computed  as: 

Phys  =~p  g  Z  (4  28) 

The  dynamic  pressure  is  computed  as: 

N 

PD=YjPoi  (4.29) 


(4.26) 

(4.27) 


The  total  orbital  velocity  vector  is  given  by: 

Kota,=fK-  (4.30) 

i= 1 

The  total  orbital  acceleration  vector  is  given  by 

=  (4.31) 

i= 1 


27 


4.11  Short-Crested  Irregular  Seaway  Kinematics 

The  kinematics  for  a  short  crested  irregular  seaway  are  computed  in  a  similar  manner  to  a 
long-crested  seaway.  A  spreading  function,  G(/3),  is  introduced  in  order  to  distribute  the  energy 
of  each  wave  frequency  over  several  directions.  Two  common  spreading  functions  are  the 
cosine  squared  and  cosine  to  fourth  functions. 

Definition  of  parameters  in  spreading  function 

[)  :  Wave  component  direction 

P0  :  Seaway  direction 

P  A  :  Spreading  angle 


(4.32) 


Lower  limit  for  beta  =  P0 - — 


'  Upper  limit  for  beta  =  P0  + 
Cosine  Squared  Spreading  Function: 

GW=  2  .  *) 

—  •cos  {fi-fio)-  — 

Pa  V  Pa 

Where :  <  f  <  fiv 

Cosine  to  the  Fourth  Spreading  Function: 

8  1  4f,  .  x) 

3  PA  r  Pa) 

Where :  p  L  <p  <  p(! 

The  spectral  density  then  becomes: 

S(co,p)=  S(co)-G(p) 

The  amplitude  for  the  wave  components  is  given  by: 

=  p-ffa, 

The  wave  elevation  at  an  arbitrary  location  is  the  computed  as: 


Cro,ai-jlilk,MXt^ 


(4.33) 

(4.34) 

(4.35) 

(4.36) 

(4.37) 


The  hydrostatic  pressure  is  computed  as: 

Phys=~p-g-z  (4.38) 

The  dynamic  pressure  is  computed  as: 
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(4.39) 


Pd=YLPdu 

i= 1  7=1  1 


The  total  orbital  velocity  vector  is  given  by: 


TV  Nfi 


^ r Total  ^ i,j 

i= 1  y=l 

The  total  orbital  acceleration  vector  is  given  by: 


TV  % 


F, 


(4.40) 


(4.41) 


4.72  Disc  ret  iza  tion  of  Spectra 

The  frequencies  of  the  discretized  spectrum  may  have  equal  or  unequal  spacing  between 
them.  In  the  case  of  unequal  spacing,  the  frequencies  are  selected  so  that  the  area  of  each  slice  of 
the  spectrum  (S(co,)-  Scot)  is  constant.  This  yields  wave  component  amplitudes  that  are  constant 
across  frequencies  for  a  given  wave  direction  /?. 

The  spectra  are  truncated  so  that  wave  components  with  infinitesimal  amplitudes  or 
extremely  high  frequencies  are  excluded,  as  they  do  not  contribute  significantly  to  the  overall 
result  and  add  computational  burden.  The  truncation  limits  are  set  so  that  significant  wave 
height  computed  from  the  area  under  the  discretized  spectrum  is  less  than  0.5%.  A  lower 
truncation  limit,  a>wm,  of  0.6 -co,,,  and  an  upper  truncation  limit,  couum,  of  4.0- com  are  appropriate  for 
all  spectra  defined  in  the  Theoretical  Spectra  section. 

4.13  Theoretical  Spectra 

There  are  several  commonly  used  theoretical  spectra  that  are  useful  in  engineering  analysis. 
The  most  widely  used  spectra  for  ship  design  evaluations  are  Bretschneider,  JONSWAP,  and 
Pierson-Moskowitz. 


4.13.1  Pierson-Moskowitz  -  Input  Hs  only 


Spectral  definition: 


Modal  period: 


S(cu) 


0.0081  -g2 


( 

•  exp 

v 


-  0.032 -g2^ 
Hs-co4  , 


(4.42) 


(4.43) 


4.13.2  Bretschneider  -  Input  Hs  and  ©m 

Spectral  definition: 
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(4.44) 


s(®)=^'4'».2'«p 

4  co 

4.1 3.3  JONSWAP  -  Input  Hs  and  com 

Spectral  definition: 


-1.25  - 


/  \4 
'  co„  ' 


V  M  ) 


where, 


JONSWAP  5 


co 


exp 


4 

-1.25 

V  co  ) 

r 


exp  -(co-(om  f  !Wo>l ) 


A  =  —  ca)n  (o. 06 5 3  3y<im 1 5  +  0. 13467)"1 
16 


_  f  0.07  for  co  <  ml 
[0.09  for  co  >  com  J 
y  =  3 .^(peaking  factor) 


(4.45) 


(4.46) 


4.13.4  Ochi  Hubble  Six  Parameter  Spectrum  and  other  Bimodal 
Spectra 

Wave  spectra  obtained  from  measured  data  in  the  ocean  often  have  double  peaks  (bimodal 
frequencies).  The  Ochi  Hubble  6-Parameter  spectral  formulation  as  described  in  Ochi  and 
Hubble  (1976)  is  actually  the  sum  of  two  three-parameter  spectra,  which  in  general  can  represent 
a  spectrum  with  two  peaks.  The  six  input  parameters  are  the  modal  frequency  com,  significant 
wave  height  Hs  and  peaking  factor  /for  each  of  the  two  spectra.  Note  that  Ochi  and  Hubble 
(1976)  use  the  symbol  X  for  the  peaking  factor  instead  of  y.  Since  A  is  used  for  wavelength  we 
have  used  y,  which  is  consistent  with  the  symbol  used  for  a  peaking  factor  in  the  JONSWAP 
spectra.  The  spectral  definition  is: 


S(eo) 


4 


f 


yj 

J 


4y,+l 

CO  J 


exp 


(  4  yj  + 1 " 

( co  f 

mj 

4 

l  4  J 

l  ®  J 

(4.47) 


where  Ais  the  Gamma  function.  The  total  significant  wave  height  for  the  Ochi  Hubble  Spectrum 
is: 

H,  =  fHl+Hl  (4.48) 

Each  of  the  two  spectra  making  up  the  total  spectrum  reduces  to  the  Bretschneider  spectrum 
when  y-=  1 ,  since  T\  1)  =  1 . 

The  Ochi  Hubble  six  parameter  spectrum  is  a  common  type  of  bimodal  spectrum.  Bimodal 
spectra  can  also  be  fonned  by  summing  any  two  single  peak  spectra;  for  example,  a  two  peak 
JONSWAP  spectrum  could  be  fonned  as  the  sum  of  two  JONSWAP  spectra  with  different 
modal  frequencies.  Bimodal  spectra  can  be  used  to  represent  a  combination  of  a  wind  driven  sea 
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and  a  swell.  Swell,  which  is  composed  of  waves  generated  from  an  old  storm  traveling  across 
the  ocean,  typically  consists  of  long-crested  waves  with  longer  wave  periods.  The  direction  of 
the  swell  may  be  different  than  that  of  the  wind  driven  sea.  For  the  case  of  a  bimodal  spectrum 
composed  of  a  wind  driven  sea  and  a  swell,  the  two  spectra  forming  the  total  spectrum  may  each 
have  a  separate  primary  wave  direction,  ffa,  as  well  as  a  different  spreading  function  G(/f). 

During  the  original  storm  large  that  creates  swell  conditions,  wind  driven  waves  are 
generated  with  the  wave  spectrum  that  probably  resembles  a  fully  developed  Pierson-Moskowitz 
spectrum,  consisting  of  many  waves  over  a  wide  period  range.  As  the  waves  travel  away  from 
the  center  of  the  storm,  the  longer  waves  will  travel  faster  than  the  shorter  waves,  and  several 
hundred  miles  from  the  location  of  the  original  stonn  the  waves  will  have  separated  into  groups 
with  different  wavelengths  (note  that  wavelength  is  directly  proportional  to  wave  period).  All  of 
the  waves  will  decrease  in  amplitude  as  they  travel  out  of  the  storm  area  into  areas  with  less 
wind,  but  the  longer  waves  will  retain  more  energy,  and  most  of  the  shorter  waves  (periods  <12 
seconds)  will  decay  completely  and  die  out  within  a  few  hundred  miles  of  the  original  storm.  The 
remaining  longer  waves  will  appear  as  a  nearly  long-crested,  uniform  wave  train,  which  is  what 
is  referred  to  as  swell.  An  observer  a  thousand  miles  or  more  from  the  original  storm  location 
will  first  see  waves  of  smaller  amplitude  with  longer  periods  of  around  20  seconds.  These  are  the 
waves  from  the  low  frequency  end  of  the  original  storm  wave  spectrum  that  travel  the  fastest  and 
reach  the  observer  first.  Over  time  the  amplitude  of  the  swell  will  increase  and  the  period  will 
decrease.  For  a  typical  swell  the  most  wave  energy  is  usually  contained  in  waves  with  a  15-17 
second  period,  corresponding  roughly  to  the  modal  period  of  the  original  storm.  As  the  period 
continues  to  decrease  past  14  seconds  the  swell  amplitude  will  likely  decrease  and  the  swell  will 
eventually  fade  away. 

4.14  Spectral  Moments 

The  nth  moment  of  a  seaway  spectrum  is  calculated  as  follows: 

oo 

mn  =  |  co”  ■  S^co)  •  dco  (4.49) 

o 

The  spectral  moments  are  useful  in  the  calculation  of  several  parameters  related  to  the  seaway. 

If  a  narrow-banded,  stationary,  Gaussian  process  with  zero  mean  is  assumed,  then  the  significant 
wave  height  is  given  by: 

H=  4-V^  (4.50) 


The  average  period  of  the  seaway  is  given  by: 


T  —  0  TT  ^ 

TAVg-  2-rc - 

/«! 


(4.51) 


The  zero-crossing  period  of  the  seaway  is  given  by: 
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(4.52) 


4.15  3rd  Order  Stokes’  Wave  -  Infinite-Depth 

The  equations  listed  in  the  previous  sections  apply  to  linear  waves.  In  this  section  the 
equations  for  a  monochromatic  2nd  and  3rd  order  Stokes’  wave  in  infinite  water  depth  will  be 
presented.  Stokes’  waves  are  based  on  nonlinear  solutions  for  plane  progressive  waves  derived 
from  power  series  expansions  in  the  wave  amplitude.  Stokes’  waves  of  varying  orders  can  be 
derived  by  retaining  more  terms  in  the  expansion. 

For  a  2nd  and  3rd  order  Stokes’  wave  in  deep  water,  the  dispersion  relation  becomes: 


co2  =gk(l  +  k2a2)  (4.53) 

where  wave  number  is  defined  as: 

k  =  ^  (4.54) 

Note  that  many  of  the  simple  relationships  between  the  wavelength  and  the  wave  frequency 
that  are  valid  for  linear  waves  (Eqns.  ((4.1)-  (4.4))  were  derived  from  the  simpler  linear 
dispersion  relation,  and  cannot  be  used  for  higher  order  Stokes’  waves.  The  relationships  shown 
in  Equations  (4.54)  through  (4.56)  are  valid  for  a  monochromatic  Stokes'  wave  of  any  order. 


Frequency: 


Wave  Celerity  (Phase  Velocity): 


(o  =  2-n  ■  f  = 


2  •  n 
T 


c  = 


CO 

k 


(4.55) 


(4.56) 


If  the  wavelength  or  wave  number  is  specified,  the  wave  frequency  and  period  can  be 
computed  from  Eqns.  ((4.53)-(4.55)).  If  the  frequency  or  period  is  specified,  Eqn.  (4.53)  must  be 
solved  to  obtain  the  wave  number.  This  equation  has  two  complex  roots  and  only  one  purely  real 
root.  The  solution  for  the  wave  number,  k,  for  a  specified  frequency  is: 

bag  ax^  (4.5 

g2+27coV)r 

When  the  nonlinear  dispersion  relation,  Eqn.  (4.53),  is  used  in  place  of  the  linear  dispersion 
relation,  the  linear  equations  for  the  velocity  potential,  Eqns.  (4.1 1)  and  (4.12),  are  accurate  up  to 
and  including  terms  of  order  ( ka )3. 

In  order  to  simplify  the  equations,  the  following  substitution  will  be  made: 


where, 


X  =  (l08ftf  a  +  12V3-J4 
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0  =  cot-  kX cos(fi)  -  kY sin(fi)  +  e 


(4.58) 


The  wave  elevation  for  a  2nd  order  Stokes'  wave  is: 


C  =  acos(@)  +  ^-ka2cos(26>) 


(4.59) 


The  wave  elevation  for  a  3rd  order  Stokes'  wave  is: 


C  =  acos(6>)  +  ^ka2cos(2@)  +  ^k2a3cos(36>) 


(4.60) 


For  linear  and  2nd  order  Stokes'  waves,  the  peak  to  trough  wave  height,  H,  is  equal  to  2A. 
For  a  3rd  order  Stokes'  wave,  the  wave  height  is  expressed  as: 


H  =  2a  H — k\k 
4 


(4.61) 


Note  that  for  anything  other  than  linear  waves,  the  amplitude,  A,  no  longer  relates  to  a 
physical  dimension,  such  as  the  distance  from  the  mean  water  level  to  the  peak  or  trough. 

Since  the  velocity  potential  used  for  the  linear  wave  is  valid  through  third  order  when  the 
wave  number  and  frequency  satisfy  the  nonlinear  dispersion  relation  in  Eqn.  (4.53),  the  wave 
orbital  velocity  can  be  obtained  by  differentiating  the  fonnula  for  the  velocity  potential.  The 
equations  in  (4.16)  use  the  linear  dispersion  relation  and  are  valid  only  for  linear  calculations. 
The  correct  equations  for  the  orbital  velocity  accurate  up  to  and  including  terms  of  order  (ka)3 
are  given  in  Equations  (4.62)  below: 


u,„  =  ■ 


v  =  {uw,vw,ww}=  vst(p)  = 
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w,„  = 


CO 


cos (co- 1  -  k  ■  X  ■  cos(7?)  -k-Y  ■  sin(/?)  +  s)  ■  sin( f5) 
ekZ  ■  sin(<y  t-k  ■  X  ■  cos  (/?)  -  k-Y  ■  sin  (j3)  +  s) 


The  pressure  at  any  point  under  the  wave  can  be  computed  from  the  Bernoulli  Equation.  In 

-  2 

deriving  the  first  order  equation  for  dynamic  pressure,  Eqn.  (4.15),  the  V  term  in  the  Bernoulli 
Equation  was  not  included  as  this  was  of  the  order  ( ka)~ .  For  2nd  and  3rd  order  Stokes'  waves 
this  term  must  be  retained.  The  total  incident  wave  pressure  for  a  2nd  and  3rd  order  Stokes' 
wave,  including  both  the  hydrostatic  and  dynamic  components  is: 


p  =  -pgZ  +  pgaekz  cos(0)  -  —  p 
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gak 
y  co  ) 


J-kZ 


(4.63) 


The  total  pressure  specified  by  the  above  equation  approaches  zero  at  the  water  surface,  so 
it  is  not  necessary  to  apply  the  pressure  stretching  used  in  the  linear  wave  theory.  Since  the 
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formulae  for  the  velocity  potential  of  a  monochromatic  2nd  and  3rd  order  Stokes'  wave  are 
identical,  and  do  not  include  terms  of  O (ka)2  or  higher,  the  expressions  for  the  wave  orbital 
velocity  and  pressure  given  in  Equations  (4.62)  and  (4.63)  are  the  same  for  both  2nd  and  3rd 
order  Stokes'  waves.  Only  the  expressions  for  wave  elevation  are  different  between  the  2nd  and 
3rd  order  Stokes'  waves. 


4.16  Inputs  and  Defaults 

Tempest  should  be  able  to  use  an  arbitrary  number  of  linear  seaways.  There  are  three  ways 
to  specify  the  input  parameters  for  a  seaway  and  each  seaway  should  be  able  to  use  a  different 
method.  The  input  methods  and  their  associated  options  are: 

•  Regular  wave  inputs 
°  Required: 

■  Wave  amplitude:  A 

■  Either  period,  length  or  frequency:  T,  X,  or  co 
0  Optional: 

■  e  -  default:  0  deg 

■  Order  of  Stokes'  wave  -  default:  linear 

•  Selection  of  a  theoretical  spectrum 

0  Required:  Spectrum  type,  Hs,  Tm,  N 
0  Optional: 

■  Po  -  default:  0  deg 

t ouim  -  default:  0.6-  com 
cOuiim  —  default:  4.0-  com 

•  r£  (random  seed  for  automatically  selected  phase  angles)  -  default:  random 

■  Pressure  stretching  flag  -  default:  use  pressure  stretching 

■  Frequency  spacing  method  -  default:  constant  area 

•  Options  are  constant  8co  or  constant  area 

■  Spreading  function  -  default:  no  spreading 

•  Spreading  function  type  -  default:  cosine  squared 

•  Pa~  default:  90  degrees 

•  Np-  default:  5  (not  including  end  points  of  spreading  function  -  they  are 
0) 

•  Input  of  spectral  ordinates 
°  Required:  o)„  S(co,) 

0  Optional: 

■  Pi  -  default:  0  deg 

■  Si  -  default:  random 

.  mutually  exclusive  with  rE 

•  re—  default:  random 

.  mutually  exclusive  with  s, 

■  Pressure  stretching  flag  -  default:  use  pressure  stretching 

•  Input  of  wave  components 
0  Required:  0)t,  a,- 

0  Optional: 
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■  [i,  -  default:  0  deg 

■  Sj  -  default:  random 

.  mutually  exclusive  with  rs 
•  re—  default:  random 

.  mutually  exclusive  with  st 

■  Pressure  stretching  flag  -  default:  use  pressure  stretching 

4.17  Intermediate  Output 

Tempest  shall  report,  in  some  manner,  Hs  computed  from  the  discretized  spectrum.  The 
graphical  Tempest  client  shall  have  the  ability  to  display  an  autocorrelation  function  for  the 
discretized  spectrum  computed  by  a  discrete  cosine  transform  of  the  spectrum. 

5  Autopilot  Control 

Tempest  Level-0  will  include  a  simple  PID  control  algorithm  to  control  the  rudder 
deflection  and  propeller  RPM.  The  rudder  control  algorithm  will  be  set  up  to  deflect  the  rudder 
in  an  unsteady  manner  to  maintain  a  constant  heading  angle  specified  by  the  user.  In  the  Level-0 
theory,  the  propeller  RPM  is  assumed  to  be  constant,  but  a  PID  control  algorithm  is  included  for 
the  propeller  to  help  compute  the  propeller  RPM  required  to  achieve  a  desired  speed  in  calm 
water. 

The  PID  controller  for  the  rudder  requires  the  following  user-supplied  inputs: 

•  Yd~  desired  heading  angle 

•  Gp  —  gain  coefficient  for  proportional  tenn 

•  Gi  -  gain  coefficient  for  integral  tenn 

•  Gd~  gain  coefficient  for  derivative  term 

•  8max  -  maximum  rudder  deflection  angle 

•  8-  rudder  slew  rate 

•  PB  —  Proportional  band  of  controller 

•  Imax  -  upper  limit  on  the  amplitude  of  the  integral  error 


The  commanded  rudder  deflection  angle  at  time  t  is  computed  as: 

8C  =80  +  Gp{yd  -  y)+  G, J0Vz>  -  ¥^))dr  +  Gdy/  ,  (5.1) 

where  xj/  is  the  time  derivative  of  the  ship  heading,  which  is  computed  using  finite  differencing 
of  the  ship  heading  angle  at  the  current  and  previous  time  step.  The  integral  gain  tenn  above  is 
limited  by  the  user-supplied  value,  Imax: 


j0Vz> 


</„ 


(5.2) 
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To  compute  the  actual  rudder  deflection  angle  at  each  time  step,  it  is  assumed  that  the  rudder 
servo  is  a  linear  first-order  lag  with  slew  rate,  8 ,  and  proportional  band,  PB.  The  actual  rudder 
deflection  is  then  computed  as: 


8 


80  + 


1-exp 


(<?.  -<?») 


(5.3) 


where  5o  is  the  rudder  angle  at  the  beginning  of  the  time  step  and  the  second  tenn  is  the 
increment  of  the  rudder  angle.  The  increment  is  limited  by  the  rudder  slew  rate: 


At 


<8  , 


(5.4) 


and  the  total  rudder  deflection  is  limited  by  the  maximum  rudder  angle,  5max: 

\8\<8 

|  |  max 


(5.5) 


6  Forces 

Tempest  requires  the  total  force  and  moment  vector  to  be  computed  at  each  time  step  in 
order  to  solve  the  equations  of  motion  for  the  ship.  In  the  Level-0  theory  the  forces  and 
moments  are  referenced  to  the  center  of  gravity  of  the  ship.  In  this  document,  the  term  “forces” 
will  hereafter  be  assumed  to  include  both  forces  and  moments.  In  Tempest  it  will  be  assumed 
that  the  total  force  can  be  decomposed  into  various  force  contributions,  and  the  total  forces  on 
the  hull  can  therefore  be  obtained  by  summing  the  various  force  contributions.  In  the  Level-0 
theory,  the  following  force  contributions  will  be  considered: 

•  Hydrostatic  forces  -  Fhys 

•  Froude-Krylov  forces  -  Ffk 

•  Radiation  forces  -  Fraj 

•  Diffraction  forces  -  F^f 

•  Hull  resistance  -  Fresist 

•  Bilge  keel  forces  -  FBk 

•  Bare  hull  maneuvering  forces  (includes  skeg  force)  -  Fman 

•  Rudder  and  fin  forces  -  FR 

•  Propulsion  forces  -  FP 

•  Weight  of  ship  -  W 

The  force  from  the  weight  of  the  ship  is  described  along  with  the  equations  of  motion  in 
Section  3.3.  The  remaining  force  components  are  discussed  in  Sections  6.1  through  6.7.  As  the 
total  force  is  computed  as  the  linear  superposition  of  component  forces,  special  care  is  taken  to 
avoid  either  missing  important  force  components  or  the  double  counting  forces.  An  example  of 
double  counting  forces  would  be  low  frequency  radiation  forces  and  maneuvering  forces. 
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Another  example  would  be  computing  and  including  the  rudder  force  when  the  hull  maneuvering 
force  is  based  on  empirical  data  that  includes  the  rudder  force  already.  The  areas  where  double 
counting  of  forces  is  a  concern  and  how  double  counting  is  avoided  is  discussed  in  the  sections 
describing  each  force  component. 

The  total  force  in  the  ship-fixed  reference  frame  acting  on  the  ship  at  the  center  of  gravity 
of  the  ship  is  the  sum  of  the  component  forces. 

F total  =  Fhys  +  Ffk  +  Frad  +  Fdy  +  Fresist  +  Fbk  +  Fman+FR+FP+W  (6.i) 

Note  that  in  the  Level-0  equations  of  motion  the  weight  of  the  ship  is  already  accounted  for 
and  the  added  mass  terms  in  the  radiation  force  are  moved  to  the  left  hand  side.  Therefore  the 
force  vector  used  to  fonn  the  right  hand  side  of  the  equations  of  motion  described  by  Equation 
(3.27)  is  computed  as: 


F  =F  +F  +F 

DUO  -1  .r-  1  J  £T  V  1  -*•  , 


RHS 


hys 


FK 


rad  _  damping 


+  F  ,-r  +  F._,  +  F DLT  +  F . .  +  F, 


dif 


BK 


+  FP+W 


(6.2) 


The  terms  forming  the  right  hand  side  of  Equation  (6.2)  are  defined  by  Equations  (6.3), 
(6.16),  (6.23),  (6.28),  (14b),  (6.72),  (6.103)  and  (6.123). 


6. 1  Hydrostatic  and  Froude-Krylov  Force 

The  hydrostatic  pressure  force  on  the  hull  is  computed  by  integrating  the  hydrostatic 
pressure  over  the  wetted  hull  surface.  The  Froude-Krylov  force  is  computed  by  integrating  the 
dynamic  wave  pressure  over  the  wetted  hull  surface.  In  both  cases  the  wetted  hull  surface  is 
defined  as  the  hull  surface  below  the  incident  wave  surface.  In  other  words,  the  diffracted  wave 
is  not  considered  in  computing  the  hydrostatic  or  Froude-Krylov  force  on  the  hull  in  the  Level-0 
theory.  The  combined  hydrostatic  Froude-Krylov  force  is  the  integration  over  the  ship’s  surface 
of  the  total  pressure  that  would  exist  in  the  incident  wave  field  in  the  absence  of  the  ship. 

In  Tempest,  the  ship  hull  geometry  will  be  described  as  a  meshed  surface  defined  by  the 
vertices  that  define  the  mesh.  The  first  step  is  to  compute  the  total  wave  height  near  the  ship  in 
order  to  determine  which  vertices  are  wet  and  which  are  dry.  The  total  wave  height  at  any 
location  at  an  instant  in  time  can  be  computed  from  Equation  (4.27).  By  computing  the  total 
wave  height  at  each  vertex  in  the  mesh,  the  intersection  of  the  incident  wave  surface  with  the 
ship  hull  can  be  determined.  The  next  step  is  to  compute  the  hydrostatic  pressure,  phys,  and  the 
dynamic  wave  pressure,/?/),  at  each  vertex  at  the  current  time  step.  The  hydrostatic  pressure  at 
any  point  is  computed  using  Eqn.  (4.28).  The  dynamic  pressure  at  any  point  at  a  given  instance 
in  time  is  computed  using  Eqn.  (4.29).  When  computing  the  dynamic  pressure,  Wheeler’s 
method  of  pressure  stretching  is  applied  as  described  in  Section  4.8  to  ensure  that  the  total 
pressure  on  the  water  surface  goes  to  zero.  The  combined  hydrostatic  and  Froude-Krylov  force 
is  then  computed  as  the  integral  of  the  hydrostatic  and  Froude-Krylov  pressures  over  the  wetted 
hull  surface. 

Ffrys+FK  =  L  (Phys  +PD)dS  (6'3) 

"  ^wet 

The  integral  shown  in  Equation  (6.3)  is  evaluated  numerically  using  the  meshed  hull 
surface  and  the  values  of  the  hydrostatic  and  dynamic  pressure  computed  at  each  vertex. 
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6.2  Hydrodynamic  Perturbation  Force 

The  total  force  acting  on  the  ship  will  be  computed  by  linear  superposition  of  force 
contributions  from  various  sources.  The  hydrodynamic  perturbation  forces  are  decomposed 
into  radiation  forces  and  wave  exciting  forces.  The  radiation  forces  represent  the  forces  due  to 
the  oscillating  motion  of  the  ship  in  calm  water.  The  wave  exciting  forces  are  further  broken 
down  into  the  Froude -Krylov  forces  and  the  diffraction  forces.  The  Froude-Krylov  forces,  which 
are  described  in  Section  6.1,  are  the  component  of  the  wave  exciting  force  resulting  from  the 
integration  over  the  ship’s  surface  of  the  incident  wave  pressure  field  that  would  exist  in  the 
absence  of  the  ship.  The  diffraction  forces  represent  the  force  due  to  the  diffraction  of  the 
incident  wave  due  to  the  presence  of  the  ship. 

6.2.1  Radiation  Force 


6.2. 1.1  Summary  of  Theory 

This  section  describes  the  theory  for  computing  the  radiation  forces  in  the  Level-0 
implementation  of  the  time  domain  dynamic  stability  code  Tempest.  The  Level-0  Tempest  is 
intended  to  be  a  code  based  on  simple  theories  that  will  provide  an  initial  framework  into  which 
more  accurate  theories  can  be  implemented.  In  the  Level-0  theory,  the  radiation  force  at  each 
time  step  will  be  computed  from  the  added  mass  and  damping  coefficients  interpolated  at  the 
instantaneous  encounter  frequency,  and  the  “memory”  force  relating  to  the  radiation  force  will  be 
ignored. 

The  total  force  acting  on  the  ship  will  be  computed  by  linear  superposition  of  the  forces 
contributed  from  various  sources.  The  seakeeping  forces  are  decomposed  into  radiation  forces 
and  wave  exciting  forces.  The  radiation  forces  represent  the  forces  acting  on  the  ship  as  it 
undergoes  prescribed  oscillatory  motions  in  calm  water.  The  radiation  forces  are  further 
decomposed  by  examining  the  oscillatory  motion  in  each  of  the  six  degrees  of  freedom 
separately.  For  steady-state  oscillatory  motion  at  a  single  frequency,  the  hydrodynamic 
radiation  force  acting  on  the  ship  can  be  expressed  as: 

Fradi  (t)  =  Re  je!a* 

where,  £.is  the  complex  amplitude  of  the  oscillation  in  the yth  mode,  and  Ty  is  the  complex 

transfer  function  which  represents  the  force  in  the  /th  direction  due  to  a  unit  amplitude  oscillation 
in  the /h  mode,  which  can  be  decomposed  into  its  real  and  imaginary  parts: 


Em 


j= i 


(6.4) 


Tij 


(o2A..  -  icoB 


(6.5) 


where  Ay  and  By  are  real  coefficients  with  Aj  being  the  added  mass  in  the  /'lh  mode  due  to  unit 
acceleration  in  the /h  direction  and  By  being  the  damping  coefficient  in  the  7th  mode  due  to  unit 
velocity  in  the  jth  direction.  The  coefficients  Ay  and  By  are  in  general  dependent  upon  the 
frequency  of  oscillation  and  the  forward  speed  of  the  ship. 


As  the  ship  undergoes  forced  oscillation  it  will  generate  waves  that  will  radiate  outwards 
on  the  free  surface.  When  the  oscillatory  motion  is  unsteady  or  has  not  reached  steady  state,  the 
presence  of  the  non-uniform  radiated  waves  on  the  free  surface  results  in  a  time  dependent 
radiation  force.  This  radiation  force  depends  on  the  previous  motion  history  of  the  ship,  as  the 
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waves  produced  by  the  ship  will  continue  to  influence  the  pressure  in  the  fluid  and  consequently 
the  force  on  the  ship,  as  they  radiate  away.  When  the  solution  is  computed  in  the  time  domain, 
this  “memory”  force  is  represented  by  the  convolution  integral  in  the  formula  below. 


Fradi(t) 
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j= i 


—A-  x  —  B  x 

ijco  j  ijco  j 


\Kij{t-T)xj{T)dT 


(6.6) 


The  convolution  in  the  above  equation  can  be  expressed  in  terms  of  either  the  velocities  i. 
or  the  accelerations 


L  l 

J K.j (t  -  r)Xj ( z)d t  =  J L.j ( t  -  r)Xj(r)d r 


(6.7) 


where  the  impulse  response  functions  used  for  the  kernels  are  related  to  each  other  by 
differentiation: 


M0  = 


(6.8) 


The  impulse  response  functions  can  be  computed  from  the  frequency  domain  added  mass 
and  damping  coefficients  by  Fourier  sine  or  cosine  transforms  as: 


2  r°°r  l  2  p^B  Xcd) 

LJt )  =  —  [  \AJcq)  - Aii(co)\cos(cot)d(d  -  —  [  — — -sin (cot)dco 

71 F  J  j[  Jo  co 

2  2  (6-9) 

Kij{t)  =  —  JoC°[s..(fi;)-S..(co)]cos(ft>0^^  =  —\o  coIa^o))-  A{cc)\m{a)t)do) 

In  the  Level-0  implementation  of  Tempest  the  “memory”  force  tenn  in  Equation  (6.6)  will 
be  ignored.  The  user  will  input  values  for  Ay  and  By  over  a  range  of  frequencies  at  zero  speed. 
The  code  will  then  interpolate  to  detennine  the  values  of  Ay  and  By  corresponding  to  the 
encounter  frequency  based  on  a  time  averaged  value  of  the  ship  speed,  the  instantaneous 
heading,  and  the  modal  frequency  of  the  incident  wave  spectrum.  Setting  the  “memory”  force  to 
zero  is  equivalent  to  assuming  that  the  added  mass  and  damping  coefficients  do  not  vary  with 
frequency.  Based  on  this  assumption,  Atj(co)  -  A( oo)  =  0  and B.  ^co)  -  B{ oo)  =  0  ,  the  memory 

force  term  in  the  expression  for  the  radiation  force  will  be  zero.  In  the  simplified  Level-0  theory, 
the  expression  for  the  radiation  force  becomes: 


Fradl  (t)  =  X 

7=1 


(6.10) 


where  Ay  and  By*  are  the  values  of  Ay(co)  and  By( co)  interpolated  at  the  encounter  frequency. 

In  the  linear  hydrodynamic  theory  used  for  Level-0,  the  added  mass  and  damping 
coefficients,  Ay  and  By,  are  computed  with  the  ship  in  an  upright  position  on  the  calm  free 
surface,  with  the  assumption  that  the  motions  are  small.  This  results  in  a  paradox  when  applied 
to  cases  where  the  ship  motions  result  in  large  pitch  and/or  roll  angles.  Consider  a  case  where 
the  ship  roll  angle  is  close  to  90°  and  the  ship  is  forced  to  oscillate  in  the  ship-fixed  z-direction. 
In  this  case  the  ship  is  “heaving”  with  respect  to  the  ship-fixed  axis  system  but  “swaying”  with 
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respect  to  the  upright  axis  system.  If  A33  and  B33  are  used  to  compute  the  force  in  ship-fixed  z- 
direction,  the  force  will  be  based  on  calculations  made  with  an  upright  ship  oscillating  vertically 
on  the  free  surface.  If  A22  and  B22  are  used,  the  force  will  be  based  on  calculations  made  with  the 
upright  ship  oscillating  horizontally  on  the  free  surface.  When  applying  linear  theory  to  a  case 
with  a  large  pitch  and/or  roll  angle,  a  choice  must  be  made  whether  to  use  coefficients  computed 
with  approximately  the  correct  geometry,  but  with  the  incorrect  motion  relative  to  the  free 
surface,  or  to  use  coefficients  computed  with  an  incorrect  geometry,  but  with  the  correct  motion 
relative  to  the  free  surface.  For  consistency  with  the  ship-fixed  mass  matrix,  the  Ajj  and  Btj 
tensors  are  transformed  from  the  upright  coordinate  system  (yawed  earth-fixed)  to  a  ship-fixed 
coordinate  system,  which  is  applied  in  the  Level-0  implementation  of  Tempest. 


6.2. 1.2  Input 

The  user  will  provide  a  set  of  added  mass  and  damping  coefficients  for  the  ship  at  zero 
speed.  Only  the  non-zero  coefficients  listed  at  the  end  of  this  section  will  be  provided,  and  these 
will  be  provided  over  a  sufficient  range  of  frequencies  to  allow  the  program  to  interpolate  values 
at  the  instantaneous  encounter  frequency.  The  transfer  functions  should  be  defined  in  the 
coordinate  system  with  x  forwards,  y  to  port  and  z  upwards  with  the  origin  at  the  ship  center  of 
gravity.  If  a  different  origin  is  used  to  define  the  transfer  functions,  the  origin  location  must  also 
be  specified  as  input  and  Tempest  will  perform  a  transformation  so  the  moments  are  referenced 
to  the  center  of  gravity.  A  wave  frequency  corresponding  to  the  modal  period  will  also  be 
specified  for  computing  the  encounter  frequency.  The  Level-0  theory  will  assume  that  the  ship  is 
symmetric  port  and  starboard.  With  this  assumption  the  motions  and  forces  in  the  lateral  plane 
and  vertical  plane  can  be  considered  separately,  and  any  term  with  one  odd  index  and  one  even 
index  will  be  zero. 
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4>12  —  45 14  -4516  -  4521  -  4523  “  4525  “  4532  “  4534  “  4536  ~  U 

n  0  —  d  0  —  d  0  —  n  0  —  n  0  —  n  0  —  n  0  —  n  0  —  n  0  —  rv 

# 41  “  # 43  ~  #45  ~  #52  “  #54  “  #56  “  #61  “  #63  “  #65  “  U 


The  “0”  superscript  is  used  to  indicate  zero  speed.  At  zero  speed,  the  added  mass  and  damping 
tensors  are  symmetric: 


4  =A 


0  B°=B°. 

i  ’  V  Jl 


(6.12) 


This  results  in  only  12  unique  non-zero  added  mass  and  damping  coefficients  that  the  user  must 
provide  as  input.  These  coefficients  can  be  obtained  from  either  a  strip  theory  code  or  from  a 
panel  method  such  as  WAMIT.  In  order  to  interpolate  values  of  each  coefficient  at  the  encounter 
frequency,  the  user  will  need  to  compute  the  values  over  a  range  of  frequencies  covering  the 
range  of  expected  encounter  frequencies.  The  coefficients  should  be  provided  as  dimensional 
values  with  the  units  as  specified  in  Table  6-1 : 
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Table  6-1  Units  for  Added  Mass  and  Damping  Coefficients 


i 

./ 

Units  for  Ay 

Units  for  Bu 

1,2,3 

1,2,3 

kg 

kg/sec 

1,2,3 

4,5,6 

kg  m 

kg  m/sec 

4,5,6 

1,2,3 

kg  m 

kg  m/sec 

4,5,6 

4,5,6 

kg  nr 

kg  m2/sec 

The  following  is  a  list  of  the  3-D,  zero-speed  added  mass  and  damping  coefficients  that  the 
user  will  provide  as  input: 

A 11°,  ^4 13°,  ^15°,  All,  Ai4,  Ai6°,  Aj3°,  A350,  A440,  A46°,  Ass°,  A(,6° 

#11  ,  #13  ,  #15  ,  #22  9  £>24  ,  £>26  ,  £>22  ,  #35  ,  #44  ,  #46  ,  #55  ,  #66 


6.2. 1.3  Implementation 

The  zero  speed  added  mass  and  damping  coefficients  will  be  provided  over  a  range  of 
frequencies  as  input  to  the  program.  At  each  time  step  the  program  will  interpolate  each  Ay  and 
Bj  coefficient  to  obtain  a  value  corresponding  to  the  encounter  frequency  based  on  the  modal 
frequency  of  the  incident  wave.  The  interpolated  values  will  then  be  corrected  to  account  for  the 
forward  speed  of  the  ship.  In  calm  water,  the  zero  frequency  value  of  Ay  should  be  used  and  By 
should  be  set  to  zero.  No  forward  speed  correction  is  applied  to  Ay  in  calm  water. 

The  Ay  and  Bj  coefficient  matrices  are  transformed  from  an  upright  coordinate  system  to  the 
ship-fixed  coordinate  system.  The  added  mass  terms  will  be  applied  to  the  left  hand  side  of  the 
equations  of  motion  as  described  in  Section  3.5.  The  force  from  the  damping  terms  will  be 
computed  and  added  as  a  component  of  the  forces  included  on  the  right  hand  side  of  the 
equations  of  motion. 

The  formulae  used  to  correct  the  added  mass  and  damping  coefficients  for  forward  speed 
were  taken  from  Principles  of  Naval  Architecture  (PNA  1989).  These  formulae  are  derived 
from  the  strip  theory  outlined  in  Salvensen,  Tuck,  and  Faltinsen  (1970).  The  formulae  in  PNA 
are  modified  from  those  in  the  original  Salvensen  et.  al.  (1970)  work,  in  that  the  original  theory 
included  transom  stern  corrections  that  have  been  omitted  and  the  original  work  did  not  include 
the  surge  degree  of  freedom,  which  has  been  added.  The  following  terms  do  not  have  forward 
speed  corrections  and  will  be  set  to  their  zero  speed  values: 

A\3  =^31=^13° 

B\3  =  B$\  =  B\3 

Al4  =  A42=  Al4 
Bl4  =  f>42=  B14 

An  =  An° 

Bn=Bn° 

An=  Aii 

Bn=  Bii 

A33—  A330 
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Bl3~  $33° 
A44=  A44 


B44-  B44 


The  remaining  terms  will  be  corrected  to  account  for  forward  speed  using  the  formulae 
listed  in  Table  6-2  below: 


Table  6-2  Forward  speed  corrections  for  added  mass  and  damping  coefficients. 


Vertical  Plane 

Lateral  Plane 

A  -  A  0  b 

A  5  _  ^15  2  ^13 

A  -  A  0  1  B 

^26  —  ^26  T  2  ^22 

C00 

0) 

e 

e 

A  —  A  0  +  B 

A  -  A  0  — 0  B 

A51  ~  A5  +  2  -“13 

62  ~  ^26  2  -°22 

0Je 

<°e 

Bl5=B15°+U0Au 

B26  =  B2(,  —  6,,(i  A  22 

B5i  =  B15°  -U0Au 

B62  =  B26  +  U0A22 

A  -  A  0  B 

^35  ^35  2  33 

A  -  A  0  +  B 

^M6  _  ^46  ^  2  n24 

A  —  A  0  1  B 

^53  —  ^35  T  2  ^33 

A  -  A  0  —  b 

^64  —  ^46  2  ^24 

CO 

CO 

e 

e 

B35  =  ^35  +  U0A33 

B46  =  B46  —  U0A24 

B53  =  B35  —  U0A33 

B64  =  B46  +  u0a24 

A  -  A  0  +  U°2  A 

^55  ~  ^55  ^  2  ^33 

A  —  A  0  1  ^0  A 
^66  ~  ^66  2  ^22 

b„=b1;+^tb„ 

Btt=Bst°+%Ba 

U0  is  the  mean  forward  speed  of  the  ship.  In  Tempest  (Level-0)  Uo  will  be  computed  as  the 
running  average  of  the  instantaneous  velocity  of  the  ship  along  the  ship-fixed  x-axis  over  an 
interval  of  time  before  the  current  time  step.  The  size  of  the  time  interval  used  for  the  averaging 
will  be  an  optional  user  input,  and  the  default  size  of  the  time  interval  will  be  15  seconds.  coe  is 
the  encounter  frequency,  which  is  computed  using  the  formula  below: 


coe  =  co - U0cos/u 

g 


(6.13) 


The  value  of  co  used  to  compute  the  encounter  frequency  will  be  the  frequency  of  the 
incident  wave  if  a  regular  wave  is  specified  and  the  modal  frequency  if  an  irregular  wave 
spectrum  is  specified.  The  angle  //  is  the  heading  angle.  It  will  be  computed  in  Tempest  at  each 
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time  step  as  the  angle  between  the  ship-fixed  x-axis  and  the  primary  direction  of  the  waves,  with 
180°  corresponding  to  head  seas  and  0°  to  following  seas.  In  calm  water  cases,  the  encounter 
frequency  should  be  set  to  zero  and  the  forward  speed  corrections  should  not  be  included.  Also 
in  calm  water  the  damping  terms,  By,  should  all  be  set  to  zero. 

The  code  will  then  transform  the  added  mass  and  damping  terms  from  an  upright  system  to 
a  ship-fixed  system.  The  transformation  matrix  for  converting  from  the  yawed  upright  system  to 
the  ship-fixed  system  is: 


cos  0  0 

sin^sin#  cos  <j) 
cos^sin#  -sin^ 


-sin  6 
sin  (f>  cos  6 
cos  (j)  cos  6 


(6.14) 


where  (j)  is  the  roll  angle  and  9  is  the  pitch  angle,  as  defined  in  Section  3.2.  In  order  to  rotate 
the  6x6  added  mass  and  damping  matrices,  each  3x3  quadrant  of  the  matrix  is  transformed 
separately  as  follows: 


1 

m 

rt 

i _ 

\l  A  3 

(6.15) 

A  *  A  ’  A  ’ 

^21  ^22  yi23 

=  [rms] 

AAA 

^21  ^22  ^23 

[TmsY 

A  *  A  ’  A  ’ 

_  31  ■/132  ^33 

A  A  A 

_  31  ^32  ^33  _ 

Note  that  prime  notation  is  used  in  the  case  of  the  added  mass  and  damping  coefficients  to 
denote  values  in  the  ship-fixed  coordinate  system.  Throughout  this  document,  prime  notation 
generally  refers  to  the  yawed  upright  frame  of  reference.  However,  the  user  supplied- 
coefficients,  Ay  and  By,  calculated  in  the  yawed  upright  frame  are  conventionally  referenced 
without  prime  notation.  This  inconsistency  in  notation  in  this  document  is  a  result  of 
maintaining  conventional  notation  for  the  initially-calculated  added  mass  and  damping 
coefficients. 

After  transforming  the  matrices  for  both  the  added  mass  and  damping  coefficients,  the 
transfonned  6x6  added  mass  matrix  will  be  added  to  the  left  hand  side  of  the  equations  of  motion 
in  place  of  the  infinite  frequency  added  mass  terms  discussed  in  Section  3.5.  The  radiation 
forces  due  to  the  damping  will  be  computed  using  the  equation  listed  below  and  added  to  forces 
composing  the  right  hand  side  of  the  equations  of  motion. 

6  ,  v 

Frad  _ damping :  (t)  =  ^ -  Btj ' (x;  -xj  (6-1 6) 

j= i 

where,  x  .  is  the  mean  velocity  in  the  jth  direction.  The  values  for  x .  will  be  computed  as  the 
running  average  of  the  velocity  vector  over  the  time  interval  before  the  current  time  step.  The 
length  of  the  time  interval  will  be  the  same  used  to  compute  Uo,  and  j,  will  be  equal  to  Uo.  At 
the  start  of  the  simulation  the  averaging  will  take  place  from  t  =  0  seconds  to  the  current  time, 
until  the  total  simulation  time  exceeds  the  length  of  the  specified  time  integral. 

If  it  is  desired  to  list  the  total  ship-fixed  radiation  force  as  an  output  at  each  time  step,  this 
can  be  computed  with  the  following  equation: 
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Frad . (t)  =  £ [-4 ’ 4  -  4 ’ (i;  -  4 )] 


(6.17) 


6.2.2  Diffraction  Force 


6.2.2. 1  Summary  of  Theory 

This  section  describes  the  theory  for  computing  the  diffraction  forces  in  the  Level-0 
implementation  of  the  time  domain  dynamic  stability  code  Tempest.  The  Level-0  Tempest  is 
intended  to  be  a  code-based  on  simple  theories  that  will  provide  an  initial  framework  into  which 
more  accurate  theories  can  be  implemented. 

The  diffraction  force  in  the  Level-0  Tempest  code  will  be  calculated  from  a  set  of  pre¬ 
computed  complex  transfer  functions.  The  transfer  functions  will  be  provided  for  a  range  of 
wave  frequencies  and  heading  angles.  The  transfer  functions  are  also  dependent  on  the  speed  of 
the  ship.  The  influence  of  forward  speed  will  be  included  either  through  the  user  supplying  a  set 
of  transfer  functions  for  several  speeds  and  interpolating  to  detennine  coefficients  at  the 
instantaneous  speed  at  each  time  step,  or  by  applying  forward  speed  corrections  to  coefficients 
supplied  by  the  user  only  at  zero  speed.  The  real  part  of  the  complex  transfer  function  represents 
the  component  of  the  force  that  is  in  phase  with  the  wave,  while  the  imaginary  part  of  the 
complex  transfer  function  represents  out  of  phase  component,  ft  is  assumed  that  the  waves  are 
defined  as  the  linear  superposition  of  long-crested  waves  traveling  in  the  direction  of  the  earth- 
fixed  X-axis.  The  incident  wave  height  at  the  center  of  gravity  of  the  ship  is  then  defined  by  the 
expression: 


1  v 

g(X,  t )  =  sin  (a>nt  -  knX  +  en ) 


(6.18) 


where  there  are  N  wave  components.  con,  k„,  sn,  Q,  are,  respectively,  the  frequency,  wave 
number,  phase  and  amplitude  of  wave  component  n.  Details  of  computing  the  incident  wave 
height  for  cases  where  the  wave  components  do  not  all  travel  along  the  earth- fixed  X-axis  are 
given  in  Section  4.11.  If  a  wave  component  is  traveling  at  an  angle  /?  to  the  earth- fixed  X-axis, 
{-k,X)  should  be  replaced  with  (-k,X cos/?  -  knYsinj3)  in  all  of  the  equations  found  in  this  section, 
where  X  and  Y  are  the  positions  of  the  center  of  gravity  of  the  ship  with  respect  to  the  earth-fixed 
origin  in  the  earth-fixed  X-direction. 

The  total  diffraction  force  in  the  ith  direction  will  be  the  sum  of  the  diffraction  force 
contributions  from  each  wave  component  defining  the  wave  system. 

Fdift (t)  =  X  (Rc[^" K > /')k  ~knX  +  sn )  +  Im[/4 (ojn , //)]-  cos (a>nt  -kX  +  sn )}  (6- 1 9) 


The  complex  diffraction  transfer  function  for  the  force  in  the  1  direction  is  expressed  as  a 
function  of  frequency  and  heading  angle: 


FiD(o)n,ju )  =  Re[f4  (con ,//)]+ ?  Im[f)D  (con ,  //)] 


(6.20) 
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In  the  linear  seakeeping  theory  used  for  Level-0,  the  diffraction  force  transfer  function 
coefficients  are  computed  with  the  ship  in  an  upright  position  on  the  mean  free  surface,  with  the 
assumption  that  the  wave  amplitudes  are  small.  The  resulting  forces  will  be  computed  in  a 
yawed  upright  coordinate  system  with  origin  at  the  center  of  gravity.  These  forces  must  be 
transformed  to  the  ship-fixed  coordinate  system  with  origin  at  the  center  of  gravity  before  they 
are  added  to  the  other  component  forces. 

6.2.2.2  Input 

The  user  will  provide  a  set  of  diffraction  force  complex  transfer  function  coefficients 
computed  for  the  ship  either  at  zero  speed  or  for  a  range  of  speeds.  The  set  of  transfer  functions 
must  cover  a  full  range  of  headings  and  wave  frequencies.  The  transfer  functions  should  provide 
forces  and  moments  about  an  origin  at  the  center  of  gravity  of  the  ship.  The  transfer  function 
coefficients  should  be  provided  as  dimensional  values  with  the  units  specified  in  Table  6-3: 

Table  6-3  Units  for  Diffraction  Transfer  Functions 


i 

Units  for  F[\ 

1,2,3 

N/m 

4,5,6 

N  m/  m 

The  transfer  functions  should  be  defined  in  the  coordinate  system  with  x  forwards,  y  to 
port,  and  z  upwards  with  the  origin  at  the  ship  center  of  gravity.  If  a  different  origin  is  used  to 
define  the  transfer  functions,  the  origin  location  will  be  specified  as  input  and  Tempest  will 
perform  a  transformation  so  the  moments  are  referenced  to  the  center  of  gravity. 

6.2.2.3  Implementation 

The  diffraction  force  coefficients  will  be  interpolated  twice  to  obtain  transfer  functions  at 
the  current  wave  heading  for  each  wave  frequency  used  to  make  up  the  incident  wave  spectrum. 
The  first  interpolation  will  be  to  obtain  a  set  of  transfer  functions  at  each  of  the  wave  frequencies 
used  by  the  incident  wave  spectrum.  The  set  of  interpolated  transfer  functions  will  be  a  function 
of  the  input  wave  heading  angles.  Since  the  wave  frequencies  of  the  component  waves  forming 
the  wave  spectrum  are  constant,  this  first  interpolation  can  be  performed  prior  to  the  start  of  the 
time  stepping  procedure,  with  the  interpolated  values  being  stored  and  used  at  each  time  step. 

The  second  interpolation  will  be  perfonned  at  each  time  step.  The  second  interpolation  will 
derive  the  transfer  function  at  the  current  wave  heading  angle  at  every  wave  frequency  used  in 
the  wave  spectrum. 

After  the  interpolation,  the  diffraction  transfer  functions  will  be  corrected  to  account  for  the 
forward  speed  of  the  ship.  The  formulae  used  to  correct  the  complex  diffraction  transfer 
functions  for  forward  speed  follow  the  formulae  used  in  SMP  and  as  described  in  Appendix  B  of 
Meyers,  et.  al.  (1981).  Forward  speed  corrections  are  only  applied  to  the  pitch  and  yaw  moment 
coefficients: 

F\ 

K 

fd 

1  3 

fd 

I4 


=  rl 


=  F 

r2 

=  f ; 


DO 


=  fa 


DO 
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Vo  is  the  mean  forward  speed  of  the  ship.  In  Tempest  Level-0  it  will  be  set  to  the  the 
running  average  of  the  instantaneous  velocity  of  the  ship  along  the  ship-fixed  x-axis.  coe  is  the 
encounter  frequency,  which  is  computed  using  the  formula  below: 


( oe=a> - U0  cos  /j 

g 


(6.22) 


The  value  of  coe  will  be  computed  separately  for  each  component  wave  frequency  in  the 
wave  spectrum.  The  angle  //  is  the  heading  angle.  It  will  be  computed  in  Tempest  at  each  time 
step  as  the  angle  between  the  ship-fixed  x-axis  and  the  primary  direction  of  the  waves,  with  180° 
corresponding  to  head  seas  and  0°  to  following  seas. 


Alternatively,  the  user  may  input  coefficients  for  both  the  radiation  and  diffraction  forces  at 
several  speeds.  Coefficients  including  forward  speed  can  be  computed  by  programs  such  as 
AQWA,  FD-Waveload,  and  Precal.  With  this  option  the  forward  speed  corrections  will  not  be 
applied,  but  an  additional  interpolation  will  be  required  to  determine  the  coefficients  applicable 
to  the  instantaneous  speed  of  the  ship  at  each  time  step. 

After  the  coefficients  are  interpolated  and  corrected  for  forward  speed,  the  diffraction  force 
can  be  computed. 


Fdifi (0  =  X {RehD K >  A)k  sin(^V  -  Kx  +  sn)+  Im[F" (oJn , /u)\n  cos (aHt  -  kX  +  en )}  (6'23) 


n= 1 


The  diffraction  forces  and  moments  will  then  be  transformed  from  the  yawed  upright 
coordinate  system  to  the  ship-fixed  coordinate  system.  The  diffraction  forces  are  then  added  to 
forces  composing  the  right  hand  side  of  the  equations  of  motion. 

6.3  Resistance  (Axial) 

6.3.1  Summary 

This  section  presents  a  method  for  computing  the  drag  force  on  a  ship  in  waves  (or  calm 
water)  to  be  implemented  in  the  time  domain  dynamic  stability  code  Tempest.  This  method  is 
appropriate  for  Theory  Level-0  through  Level-II.  The  Level-Ill  theory  will  capture  certain 
elements  (wave  making  drag  in  particular)  of  the  method  presented  herein.  The  last  sub-section 
of  will  cover  necessary  changes  when  the  Level-Ill  Theory  is  implemented 

6.3.2  Input 

The  user  shall  supply  a  curve  of  total  resistance  for  the  ship  defined  by  several  points.  Let 
Rt  be  the  total  resistance  of  the  ship  and  V  be  ship  speed.  The  total  resistance  coefficient,  CY,  is 
computed  according  to  the  following  equation: 
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(6.24) 


C  =  T 
T  0.5  pScwV2 

where: 

Ct  =  total  resistance  coefficient 
Rt  =  total  ship  resistance 
V=  ship  speed  in  m/s 
p  =  water  density 

Sew  =  wetted  surface  of  the  ship  in  calm  water 


6.3.3  Scaling 

If  the  simulation  is  to  be  performed  at  a  scale  other  than  that  for  which  the  resistance  curve 
is  given,  then  the  frictional  resistance  coefficient  of  the  ship  is  computed  at  each  input  speed 
according  the  ITTC  1957  guideline  as: 


Q  = 


0.075 

(Log10(Rn)-2)2 


(6.25) 


Where: 

Cf=  the  frictional  resistance  coefficient 

Rn  =  the  Reynolds  number  of  the  ship  at  the  input  speeds 

The  residuary  resistance  coefficient  (includes  all  resistance  components  other  than 
frictional)  is  then  calculated  for  each  input  speed  as: 


c  =C  -C 


(6.26) 


A  new  list  of  C/ values  is  computed  using  Reynolds  numbers  based  on  the  scaled  ship 
length,  rather  than  the  input  length.  A  new  table  of  Ct  values  is  then  constructed  using  the  scaled 
values  of  C/and  the  Cr  values  from  Equation  (6.26). 


6.3.4  Time  Domain  Calculations 

At  each  time  step,  Ct  is  interpolated.  The  velocity  used  in  the  interpolation  is  given  by: 


V  =U 

v  INT  ^  G 


(6.27) 


where: 

V iNT  =  the  speed  to  use  for  the  interpolation  of  Ct 

Ug  =  the  axial  velocity  of  the  ship  (in  the  ship-fixed  x  direction) 

The  resistance  at  each  time  step  is  then  computed  by 

F resist  =  — 0 .5C \  _int  PS jf INT  (6.28) 

where: 

F resist =  the  resistance  force  computed  for  the  time  step 

Ctjnt=  the  total  resistance  coefficient  interpolated  using  VINT  from  Equation  (6.27) 
p  =  water  density 

Si  =  the  instantaneous  wetted  surface  obtained  from  an  integration  of  panel  areas 
Vint=  the  instantaneous  speed  as  defined  in  Equation  (6.27) 
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6.3.5  Reference  Frame 

The  resistance  (Fresist)  acts  in  the  longitudinal  direction  of  the  ship-fixed  reference  frame 
and  is  applied  at  the  origin  of  the  ship-fixed  frame,  so  that  no  pitch  moment  is  induced  directly 
by  the  resistance  force.  The  sinkage  force  and  trim  moment  acting  on  the  ship  in  calm  water  due 
to  its  forward  speed  can  be  added  in  separately,  but  will  not  be  computed  using  the  resistance 
force  vector.  The  steady  sinkage  force  and  trim  moment  will  not  be  included  in  the  Level-0 
implementation  but  may  be  introduced  in  Level-11  or  Level-111.  In  that  case  the  steady  sinkage 
force  and  trim  moment  will  be  determined  by  interpolation  of  a  set  of  user-supplied  force  and 
moment  coefficients  similar  to  the  coefficients  Ct  described  in  the  input  section.  The  variable 
Vint  will  be  used  for  the  interpolation  of  these  coefficients  as  well,  but  scale  effects  will  be 
ignored  for  the  sinkage  force  and  trim  moment  coefficients. 


6.3.6  Extension  to  Level-Ill  Theory  Framework 

For  the  case  where  the  wave  making  drag  is  captured  by  the  slender  body  theory,  the 
proposed  method  could  be  used  with  some  minor  alterations.  The  user  could  still  input  the  total 
resistance  curve.  Some  preliminary  simulation  runs  are  made  at  constant  forward  speeds  and  the 
wave  resistance  ( Rw )  is  computed  followed  by  CV.  A  residuary  resistance  coefficient  (Cr) 
excluding  wave  making  drag  is  then  computed  as 

Cr  =  CT-Cf-  Cw  (6.29) 

Cf  is  then  treated  as  described  in  the  previous  section  of  this  paper  and  Cr  is  assumed 
constant  with  scale.  A  resistance  coefficient  that  shall  be  called  the  non-wave-making  resistance 
coefficient  ( Cnw )  is  computed  as: 


Cnw  -  Cr  +  Cf 

And  the  resistance  force  is  computed  as 

=-0.5  CmpSIVf 


r  int 


(6.30) 

(6.31) 


6.4  Bilge  Keel  forces 
6.4.1  Preface 

The  Level-0  Tempest  theory  will  use  the  Ikeda-Himeno-Tanaka  (IHT)  method  for 
computing  the  roll  moment  from  bilge  keel  damping  (see  Himeno  (1981)).  The  Level-0  Tempest 
will  not  include  all  of  the  components  of  roll  damping  described  in  Himeno  (1981),  but  will  only 
include  the  normal  force  bilge  keel  damping  force  and  the  hull  pressure  bilge  keel  damping 
force.  The  wave  making  damping  from  the  bilge  keels  is  not  included  in  this  model.  Wave 
making  damping  from  the  bilge  keels  is  typically  very  small  and  can  be  ignored  in  the  Level-0 
Theory.  However,  if  the  bilge  keel  geometry  is  included  in  the  added  mass  and  damping 
coefficient  calculations  used  for  the  radiation  force  model,  the  wave  making  damping  from  the 


48 


bilge  keels  will  be  captured  with  no  double-counting  of  the  damping  terms  computed  in  this 
section. 

The  description  of  the  method  in  Sections  6.4.2  through  6.4.6  is  derived  directly  from  a 
MARIN  document  describing  work  they  performed  for  the  US  Navy  to  develop  a  time  domain 
implementation  of  the  IHT  roll  damping.  The  original  IHT  method  was  developed  for  frequency 
domain  methods.  There  are  several  terms  used  in  their  discussion  that  are  not  defined  in  their 
documentation.  These  will  be  defined  here,  and  a  brief  description  will  also  be  given  of  the  user 
supplied  input  needed  to  implement  the  model.  This  document  describes  only  the  bilge  keel 
damping  portion  of  the  IHT  roll  damping  method.  There  are  also  components  of  roll  damping 
from  the  hull  frictional,  eddy  making  and  lift  damping.  Only  the  bilge  keel  damping  will  be 
included  in  the  Level-0  theory  for  Tempest. 

Definition  of  terms  not  defined  in  the  MARIN  documentation: 

Hq(x)  is  the  ratio  of  the  local  half  beam  to  the  local  draft,  Ho(x)  =  B(x)/(2  T(x)) 

B(x)  in  the  definition  above  is  the  local  Beam  at  x 

the  sectional  area  coefficient  a  is  the  sectional  area  /  (B(x)-T(x)) 


User  supplied  input 


The  following  terms  will  be  specified  at  a  series  of  stations  along  the  length  of  the  bilge 
keel  (including  the  forward  and  aft  ends  of  the  bilge  keels): 


B(x)  -  local  waterline  beam  at  x  in  meters 

T(x)  -  local  sectional  draft  atx  in  meters 

S(x)  -  sectional  area  of  the  hull  section  at  x  in  square  meters 

bBK(x)  -  bilge  keel  height  at  x  in  meters 

i~bk(x)  -  distance  from  the  center  of  gravity  to  the  bilge  keel  in  meters  at  x 


From  the  input  variables  listed  above,  the  following  two  terms  can  be  computed  at  each 
position  along  the  bilge  keel: 


H  oW 


B{x) 
2T(x)  ’ 


a(x)  = 


S(x) 

B{x)T{x) 


In  addition  to  the  input  variables  specified  at  positions  along  the  length  of  the  bilge  keel, 
the  following  global  input  variables  will  be  supplied. 

con  -  natural  roll  frequency  in  radians/second. 

LBk  ~  overall  length  of  the  bilge  keel 

In  equations  (6.50)  and  (6.54),  where  the  roll  amplitude  (AJ)  is  required,  a  value  of  10 
degrees  should  be  used. 
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The  symbols  listed  in  the  table  in  the  Nomenclature  section  apply  to  both  the  FDS  and  IHT 
models.  There  are  several  symbols  listed  in  that  table  that  are  needed  only  for  the  FDS  model. 
The  description  of  the  FDS  model  has  been  removed  from  this  document,  as  it  is  based  on 
proprietary  MARIN  experimental  data. 


6.4.2  Introduction 

The  (viscous)  bilge  keel  damping  is  defined  as  the  increment  of  roll  damping  when  bilge 
keels  are  installed.  It  includes  not  only  the  roll  damping  of  the  bilge  keels  themselves  but  also  all 
the  interaction  effects  among  the  bilge  keels,  the  hull  and  the  waves.  In  the  IHT-fonnulation,  the 
total  bilge  keel  damping  is  defined  as  the  sum  of  the  normal  force  bilge  keel  damping  and  the 
hull  pressure  bilge  keel  damping: 


M  =M  +M 

BK  lv±  BK,N  T  1V1  BK,H 


(6.32) 


6.4.3  IHT  -  Normal  force  bilge  keel  damping 

The  nonnal  force  bilge  keel  damping  coefficient  is  first  calculated  for  each  section  and  then 
integrated  along  the  ship  length,  i.e.  from  the  bilge  keel  aft  end  to  the  bilge  keel  fore  end: 

B  =  [  B  (x)-dx  (6.33) 

nBK,N  J  °BK,N\X)  UX 

xBK,aft 


The  sectional  normal  force  bilge  keel  damping  coefficient  consists  of  a  zero  speed  part  and 
a  forward  speed  part: 


BBK,N  (•*•)  BBK,N, 0  (•*•)  +  BBK,N,U  (•*•) 


(6.34) 


The  zero  speed  sectional  normal  force  bilge  keel  damping 

8  ,  ,  ,  (  22.5  (x) •  A, 

Bbk,n,0 ( X. )  =  - — ; -p-rBK O)  •  bBK O)  •  /  0) •  •  an  ■  I - + 


A 


3 -n 


x-f(x)  bBK(x)  ) 


(6.35) 


where  / (x)  is  defined  and  elaborated  on  in  the  appendix  (for  layout  reasons). 

The  forward  speed  sectional  normal  force  bilge  keel  damping  coefficient  is  expressed  as 

(6.36) 


B 


BK,N,U 


0)  =  y-p- b\K  (x)  •  r]K  (x)  •  U 


From  the  above  expressions,  we  observe  that  part  of  the  zero  speed  sectional  normal  force 
bilge  keel  damping  coefficient,  Bbk  JV0(x)  ,  is  proportional  to  the  roll  velocity;  the  other  part 

depends  on  the  roll  amplitude  only;  the  forward  speed  sectional  normal  force  bilge  keel  damping 
coefficient,  Bbk  n  u(x)  ,  is  proportional  to  the  forward  speed  and  does  not  depend  on  the  roll 

velocity. 

These  observations  allow  us  to  write  in  a  more  direct  manner 
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B BK,N  (U ,  P )  aBK,N,0  '  P  +  0BK,Nfi  (^)  +  aBK,N,U  '  ^ 

where  the  coefficients  are  given  by 


a 


BK,N,  0 


=  2.4-- 


3 -k 


■ P ■ 


3 a.  jore 

J  VW' 


rBK(X)-f  (X)'dx 


XBK  ,aft 


(6.37) 

(6.38) 


A 


BK,N,U 


22.5 

n 


8 

3-;r 


•P- 


xBK,fore 

j  b2BK<,x)-r2BK(x)- f(x)-dx 

xBK,aft 


(6.39) 


XBK,fore 

® BK,N,U  ’  P  ’  J  ^ BK  C^)  '  rBK  (•*)  ' 

XBK  , aft 


(6.40) 


6.4.4  IHT  -  Hull  pressure  bilge  keel  damping 

The  hull  pressure  bilge  keel  damping  coefficient  is  first  calculated  for  each  section  and  then 
integrated  along  the  ship  length,  i.e.  from  the  bilge  keel  aft  end  to  the  bilge  keel  fore  end: 

XBK  .fore 

B  =  g  (x)-dx  (6-41) 

u BK,H  J  l,BK,N\a' )  UA 

xBK,aft 


The  sectional  hull  pressure  bilge  keel  damping  coefficient  is  expressed  as 

bbk,h  0)  =  ~r~  •  P  •  rBK  O)  •  T2  O)  • 1 U)  •  / 2  (x)  ■  con  •  4 

3 -n 


(6.42) 


From  this  expression  we  observe  that  Bbk  h  0(x)  is  proportional  to  the  roll  velocity; 
therefore,  we  may  write  in  a  more  direct  manner 

BBK,H  —  aBK,H  '  P 


(6.43) 


where  the  linear  coefficient  aBK  H  depends  on  the  bilge  keel  parameters  only: 

4 


a 


BK,H  t  '  P  ' 
3 -7t 


j  rBK  (x)  ■  T2  (x)  • I (x)-  f2(x)  ■  dx 


(6.44) 


xBK,aft 
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6.4.5  IHT  -  Bilge  keel  damping  -  Summary  (new  approach) 

In  the  new  approach  the  bilge  keel  damping  is  calculated  as 

M  bk(U  ,  p)  =  —((%BK,N,0  '  \p\  +  PbK,N,Q  +  aBK,N,U  'U  +  Cl  BK H  •  '  P  (6.45) 


where  (6.38),  (6.39),  (6.40)  and  (6.44)  are  used.  The  bilge  keel  damping  force  includes  only  a 
roll  moment  tenn.  Therefore,  the  bilge  keel  force  vector  can  be  expressed  as: 


Fbk  =(0,0,0,  Mbk,  0,0) 

(6.46) 

6.4.6  Calculation  of  underlying  terms  in  IHT-method 

f{x )  =  1  +  0.3  -exp(-160-(l-cr(x))) 

(6.47) 

/  (x)  =  -Ax  (x)  •  C;  (x)  +  A2  (x)  •  C;  (x) 

(6.48) 

c;(x)  =  i.2 

(6.49) 

C(x)  =  22.5-  bBK^x)  1.2 

P  71 '  rBK  (v)  •  fix)  •  A^ 

(6.50) 

cD(x)  =  c;(x)-c;(x) 

(6.51) 

Ax  (x)  =  (m3  (x)  +  m4  (x))  •  m%  (x)  -mf  (x) 

(6.52) 

iix)-  ™A{X) 

3  •  (H0  (x)  -  0.2 15  •  mx  (x)) 

+  (  ‘  ^  3  2  +  mx  (x)  •  (m3  (x)  •  m5  (x)  +  m4  (x)  •  m6  (x)) 

6(1  -  0.215  •  m]  (x) ) 

(6.53) 

s0  (*)  =  0.3  •  n  ■  f{x)  ■  rBK  (x)  •  A4  + 1 .95  •  bBK  (x) 


(6.54) 


m ,  (x)  = 


R(x) 

T(x) 


(6.55) 
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(6.56) 


m2(x) 


OG 

T\x) 


m3  (x)  =  1  -  m]  (x)  -  m2  (x) 


m4  (x)  =  H0  (x)  -  mt  (x) 


m6(x)  = 


0.4 14  •  H0  (x)  +  0.065 1  •  mx  (x)  -  (0.382  +  0.0 106  •  H0  (x))  •  mx  (x) 
( H0 (x)  -  0.215  •  mfx)\  (l  -0.215  •  /n,(x)) 

s0(x)  n 


m1 (x)  = 


- m,  (x)  for  sn  (x)  >  —  •  R(x ) 

T(x)  4  1  0  4 


0 


n 


for  s0(x)<  —  -R(x) 


(6.57) 

(6.58) 

(6.59) 

(6.60) 

(6.61) 


mfx)  =  { 


m7(x)  +  0.414  ■  mx(x) 


tc 


for  s0(x)>--R(x) 


m 


■  (x)  +  fl  ■ 


1  -  cos 


'loW" 


R(x) 


n 


JJ 


mfx)  for  sg{x)<--R{x) 


(6.62) 


The  bilge  keel  radius  R(x)  and  the  distance  rBK  (x)  from  the  roll  axis  to  the  bilge  keel  are 
calculated  as  follows: 


R\x)  = 

2-T(x)'J  4 

a(x)) 

for  7/0(x)<l 

and 

T(x) 

R(x)  =  < 

T(x) 

for  H0(x)>  1 

and 

R'^>i 

T(x) 

R\x) 

otherwise 

(6.63) 


(6.64) 


Gk(x)  =  t(x)' 


|t/0(x)- 

(  1  r \ 
1  —  V2 

1  R(*) I'  ,  fj  °G 

(  1  rf 

1  —  V2 

1  R(x)  Y 

l  2  J 

1  T(x)  J  [  T  1 

l  2  J 

1  T(x)  J 

—i  1/ 2 


(6.65) 
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6.5  Bare  Hull  Maneuvering  forces 
6.5.1  Summary  of  Theory 

The  maneuvering  module  in  the  Level-0  Tempest  code  will  follow  the  ship  specific 
maneuvering  model  method  that  is  based  on  an  Abkowitz-type,  coefficient-based,  maneuvering 
model  (Abkowitz,  1962).  The  method  assumes  that  data  is  available  from  rotating  arm  tests  or 
PMM  tests.  The  empirical  coefficients  used  by  the  method  are  computed  using  the  MANSIM 
program  described  by  Kopp  (2000)  and  Kopp  (2007).  An  example  of  this  maneuvering  model 
for  a  containership  is  described  in  Son  and  Nomoto  (1981,  1982)  which  includes  all  of  the 
empirical  coefficients  needed  to  implement  the  model  for  that  ship. 

In  many  cases  model  test  data  from  either  rotating  arm  tests  or  PMM  tests  include  rudder 
tenns,  while  Tempest  has  a  separate  rudder  model  to  compute  the  forces  on  the  rudder.  In  these 
cases,  the  rudder  effects  are  embedded  in  all  of  the  measured  forces  and  moments,  not  just  the 
specific  rudder  angle,  5,  terms.  Therefore,  the  model  test  rudder  terms  are  included  in  the  hull 
forces  calculations.  To  eliminate  double  counting  of  the  rudder  force,  the  entire  calm  water 
“rudder”  force  is  calculated  using  the  procedure  outlined  in  Section  0  (with  the  wave  orbital 
velocities  set  to  zero)  where  these  forces  are  then  subtracted  from  the  hull  maneuvering  forces. 
Similarly,  if  the  empirical  data  is  based  on  model  tests  that  included  a  propeller,  the  propeller 
side  forces  will  be  embedded  in  the  measured  forces  and  moments.  If  this  is  the  case,  the  calm 
water  propeller  side  force  should  be  subtracted  from  the  hull  maneuvering  forces  using  the  same 
procedure  as  was  used  for  the  rudder  forces:  computing  the  propeller  side  force  using  the 
Tempest  propeller  model  with  the  wave  orbital  velocities  set  to  zero  and  then  subtracting  this 
from  the  hull  maneuvering  force.  The  model  test  data  may  also  include  the  skeg  forces  in  all  the 
measured  forces.  In  the  Level-0  theory,  the  skeg  force  will  be  left  in  the  maneuvering  force  and 
no  separate  skeg  force  model  will  be  applied.  In  some  instances  the  model  tests  are  performed 
with  a  “bare  hull”  model  and  therefore  do  not  include  rudder  forces  or  propeller  side  forces.  In 
these  instances  the  rudder  and  propeller  side  forces  should  not  be  subtracted  from  the  hull 
maneuvering  forces,  as  no  double  counting  of  the  rudder  force  is  present.  The  input  for  the 
maneuvering  model  must  include  a  flag  to  specify  whether  or  not  the  rudder  and  propeller  side 
forces  are  embedded  in  the  hull  maneuvering  force  coefficients. 

This  method  includes  the  option  to  include  a  wide  range  of  maneuvering  coefficients.  In 
most  cases  only  a  handful  of  the  coefficients  will  be  used.  Depending  on  how  the  regression 
analysis  is  performed,  different  coefficients  will  be  used  or  not  used  in  the  maneuvering  model. 
The  Son  and  Nomoto  containership  maneuvering  model  for  instance  does  not  include  any  odd 
quadratic  terms  (i.e.  coefficients  of  v|v|  ),  so  these  coefficients  will  all  be  set  to  zero  for  the 
containership.  Son  and  Nomoto  use  only  even  quadratic  terms  and  the  odd  terms  are  either 
linear  or  cubic.  Other  maneuvering  models  may  include  both  odd  quadratic  and  cubic  terms  in 
the  equations.  Several  terms  included  in  the  Son  and  Nomoto  maneuvering  equations  will  not  be 
implemented  as  they  are  already  computed  elsewhere  within  Tempest.  The  tenns  that  will  be  left 
out  are  the  hull  resistance  and  propeller  thrust  terms  in  the  surge  equation,  as  well  as  the  rudder 
tenns  in  all  of  the  equations.  The  Son  and  Nomoto  model  appears  to  compute  the  forces  on  the 
rudder  in  a  manner  similar  to  the  Tempest  rudder  model,  so  these  terms  are  simply  left  out. 

Since  the  rudder  terms  are  not  embedded  in  the  coefficients,  there  is  no  need  to  subtract  out  the 
Tempest  calm  water  rudder  force  from  the  maneuvering  forces  when  the  Son  and  Nomoto 
maneuvering  formulation  is  used.  Son  and  Nomoto  also  include  terms  proportional  to  roll 
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velocity  in  the  roll,  sway,  and  yaw  equations.  It  is  assumed  that  these  terms  will  be  included 
within  the  roll  damping  model,  so  they  are  removed  from  the  maneuvering  equations. 


6.5.2  Input 

The  user  will  supply  the  set  of  coefficients  used  to  compute  the  maneuvering  forces.  Each 
coefficient  will  be  specified  at  several  Froude  numbers.  At  each  time  step  the  program  will 
compute  the  maneuvering  forces  corresponding  to  the  Froude  numbers  at  which  the  coefficients 
were  specified.  The  code  will  then  interpolate  to  obtain  forces  and  moments  for  the  Froude 
number  based  on  the  instantaneous  speed  of  the  ship  at  each  time  step.  For  the  Son  and  Nomoto 
(1981)  containership  maneuvering  model  the  coefficients  are  independent  of  Froude  number,  so 
no  interpolation  is  required. 

The  list  of  coefficients  which  may  be  used  in  the  equations  for  the  surge,  sway,  roll,  and 
yaw  equations  are  listed  in  Table  6-4.  Typically  only  ten  to  twelve  terms  in  each  column  will  be 
used  for  a  given  model  and  the  rest  will  be  set  to  zero. 

For  the  heave  and  pitch  equations,  a  strip-wise  cross-flow  drag  model  is  used  which 
requires  the  following  input  at  set  of  stations  along  the  hull: 

Cff(x),  B(x),  x  specified  for  nx  longitudinal  positions  along  the  hull 

Cff  is  a  2-D  cross  flow  drag  coefficient  used  to  compute  the  viscous  damping  from  a  2-D 
hull  section  as  a  result  of  the  heave  velocity  of  the  section.  The  default  value  will  be  1.0.  B(x)  is 
the  local  beam  of  the  ship  on  the  waterline.  Both  Cff(x)  and  B(x)  will  be  specified  at  a  series  of 
longitudinal  positions  along  the  length  of  the  hull.  The  drag  on  the  2-D  section  with  a  constant 
heave  velocity  is  estimated  as  Cff(x)-0.5-p-  (w(x))2 ■  B(x)  ,  where  x  is  specified  in  the  ship-fixed 
coordinate  system  with  the  origin  at  the  center  of  gravity. 

The  coefficients  for  the  surge,  sway,  roll,  and  yaw  equations  are  derived  from  the  rotating 
ann  test  data.  The  rotating  arm  test  data  covers  a  range  of  parameters  up  to  a  moderate  heel 
angle,  yaw  rate,  and  sway  velocity.  Extrapolation  outside  of  this  range  can  lead  to  errors  in  the 
computed  maneuvering  forces.  The  user  will  supply  maximum  values  of  heel  angle,  sway 
velocity  and  yaw  rate  to  be  used  in  the  maneuvering  equations:  <j)max,  vmax,  rmax.  The  absolute 
value  of  the  roll  angle,  sway  velocity,  and  yaw  will  be  capped  at  the  user  specified  maximum 
values.  In  most  cases  the  rotating  ann  data  will  reference  the  forces  and  moments  to  the  ship 
center  of  gravity.  The  input  must  specify  the  origin  used  for  the  empirical  data.  The  coefficients 
listed  in  Table  6-4  should  be  specified  in  a  coordinate  system  with  x  forwards,  y  to  port,  and  z 
upwards.  If  the  origin  used  to  obtain  the  empirical  coefficients  is  different  from  the  center  of 
gravity  of  the  ship  in  the  Tempest  simulation,  Tempest  will  first  compute  u,  v,  etc  at  the  origin 
used  for  the  empirical  data.  Then  it  will  compute  forces  and  moments  at  the  origin  used  for  the 
empirical  data.  These  forces  and  moments  will  then  be  transfonned  to  the  Tempest  origin,  which 
in  Fevel-0  is  located  at  the  center  of  gravity  of  the  ship. 

6.5.3  Implementation 

The  following  terms  are  computed  or  defined  outside  of  the  maneuvering  force  module  at 
each  time  step.  They  are  used  as  input  to  the  maneuvering  force  calculations: 

q  Angular  velocity  about  the  ship-fixed  y-axis  in  radians 

r  Angular  velocities  about  the  ship-fixed  z-axis  in  radians 
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u,v,w  Velocity  of  the  ship  at  the  ship-fixed  reference  frame  origin  in  the  ship-fixed  x,  y 
and  z  directions  respectively. 

8  The  rudder  angle  in  radians 

p  Water  density 

(/>  Roll  angle  in  radians 

The  maneuvering  equations  for  surge,  sway,  roll  and  yaw  are  all  written  in  non- 
dimensional  form.  The  forces,  moments,  velocities  and  rates  that  appear  in  the  equations  are  all 
in  non-dimensional  form.  The  forces  and  moments  are  non-dimensionalized  using  length- 
squared  or  -cubed.  To  account  for  the  change  in  the  draft  of  the  ship  as  it  maneuvers  in  waves, 
the  ratio  of  the  instantaneous  midship  draft  to  calm  water  midship  draft  is  incorporated. 

forces:  0 .5- p-L2 -U2 -Tm\ 

3  2 

moments:  0.5 -p-L-  U  -Tm,  where 
p  =  density  of  water  kg/m 
L  =  length  between  perpendiculars 

U  =  ship  velocity,  Vu2  +  v2 

Tm  =  ratio  of  the  instantaneous  draft  over  the  still  water  draft  measured  at 
midships.  Tm  is  an  attempt  to  scale  the  forces  and  moments  by  the  submerged 
surface  area. 

The  maneuvering  equations  for  surge,  sway,  roll,  and  yaw  are  based  on  rotating  ann  test 
data  which  covers  a  limited  range  of  roll  angles,  sway  velocities,  and  yaw  rates.  For  values 
above  the  maximum  values  specified  by  the  user,  the  maximum  values  from  the  test  matrix  will 
be  used  in  the  maneuvering  equations. 

The  non-dimensional  variables  for  the  sway  velocity  and  yaw  rate  are  defined  as: 
v'=  v/U  (absolute  value  not  to  exceed  vmax/U) 

r'  =  rL/U  (absolute  value  not  to  exceed  rmaxL/U) 

The  roll  angle,  (/),  and  the  rudder  angle,  8,  are  dimensional  and  specified  in  radians.  If  the 
absolute  value  of  the  roll  angle  exceeds  (f>max,  a  roll  angle  with  magnitude  (f)max  and  the  proper  sign 
is  used  in  the  maneuvering  equations.  Note  that  here  the  prime  variables  v'  and  r\  are  non- 
dimensional,  whereas  the  variables,  v  and  r,  are  dimensional.  This  differs  from  the  convention 
used  in  the  other  sections  of  the  theory  manual,  which  use  only  dimensional  velocities. 

The  equations  used  in  Tempest  to  compute  the  maneuvering  contribution  to  the  non- 
dimensional  surge  and  sway  force  and  the  non-dimensional  roll  and  yaw  moments  are 
formulated  using  the  terms  listed  in  Table  6-4.  For  each  equation,  each  coefficient  is  multiplied 
by  the  values  listed  in  the  last  column  of  Table  6-4,  and  then  the  sum  of  all  the  terms  gives  the 
total  maneuvering  force  or  moment.  For  a  given  equation,  many  of  the  coefficients  will  be  zero. 
The  coefficients  that  are  specified  as  non-zero  will  depend  on  the  maneuvering  model  specified. 
Certain  tenns  may  always  be  zero  for  a  symmetric  ship  (such  as  odd  terms  of  v  and  r  for  the 
surge  equation),  but  these  terms  are  still  included  in  Table  6-4  for  completeness. 
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Table  6-4  Coefficients  used  by  the  Maneuvering  Force  Model  for  the  Surge  and  Sway  Force  and 

the  Roll  and  Yaw  Moment. 
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(6.66) 


Fx  =  Surlnt  +  r'-SurR  +  r\r\  SurRAR  +  (r()2-SurR2  +  ... 

Fy  =  Swylnt  +  r'SwyR  +  r'\r'\  -SwyRAR  +  (r()2-SwyR2  +  ... 

Mx'=  Rollnt  +  r'-RolR  +  r'\r'\  RolRAR  +  (r()2-RolR2  +  ... 

Mz'  =  Yawlnt  +  r'-YawR  +  r'\r'\  -YawRAR+  (r()2-YawR2  +  ... 

The  forces  will  then  be  converted  to  dimensional  form  (in  N  and  N-m)  and  added  to  the  other 
force  components. 

Fx  =  Fx'-0.5-pL2-U2-Tm  ,  Fy  =  Fy'-0.5-pL2-U2-Tm  (6  67) 

Mx  =  Mx-  0.5-p-L3- U2-  Tm  ,  Mz  =  Mz'-  0.5- p-L3-  U2- Tm 

Fx,  Fy,  Mx  and  Mz  are  already  defined  in  the  ship-fixed  reference  frame  and  are  added  to  the 
other  force  components  to  compute  the  total  force  and  moment  acting  on  the  ship  at  each  time 
step. 

The  heave  and  pitch  maneuvering  forces  are  computed  by  computing  the  cross  flow  drag  at 
a  set  of  stations  along  the  length  of  the  hull,  and  then  integrating  over  the  length  of  the  hull  to 
compute  the  heave  force  and  pitch  moment.  The  local  heave  velocity  of  a  section  located  at 
longitudinal  position  x  from  the  center  of  gravity  can  be  computed  from  the  heave  velocity  at  the 
center  of  gravity,  w,  and  the  rotational  velocity  about  the  ship-fixed  y-axis,  q: 

w<x>  =  w-‘l'x  (6.68) 

The  vertical  cross  flow  drag  force  on  the  section  is  then  computed  as: 

Drag  =  ■  Cff(x)  ■  B(x)  ■  (w(x))2  (6.69) 


The  heave  force,  Fz,  and  pitch  moment,  My,  are  then  computed  by  integrating  the  drag  over  the 
length  of  the  hull: 


Fz  =\p\ Cff (*)  ’  50)  •  {w{x)fdx  (6.70) 

My  =^p\x'Cff(x)'  B(x) '  (w(*))2  dx  (6.71) 


The  integration  is  performed  over  the  wetted  length  of  the  hull.  A  trapezoidal  rule 
integration  scheme  can  be  used.  The  heave  force  and  pitch  moment  computed  using  the  above 
formulae  are  in  dimensional  form  and  are  defined  in  the  ship-fixed  frame.  Fz  and  My  are  added 
to  the  other  force  components  to  compute  the  total  force  and  moment  action  on  the  ship  at  each 
time  step.  The  complete  bare  hull  maneuvering  force  vector  is  then  defined  as: 

Fman  =  (Fv  >  Fy  >  F=  >  Mx .  My  .  Mz  )  (6.72) 


6.5.4  Zero  Speed  Performance 

Due  to  the  fact  that  these  terms  are  divided  by  ship  speed,  non-dimensional  sway  velocity, 
v',  and  yaw  rate,  r\  have  the  potential  to  be  very  large  as  ship  speed  approaches  zero.  This  is 
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avoided  by  treating  some  terms  differently  as  ship  speed  approaches  zero.  The  uncertainty  of  the 
maneuvering  model  coefficients  at  these  speeds  warrants  this  pragmatic  approach. 

Speeds  less  than  0.005m/s  are  considered  zero  speed.  In  this  case,  Fn,  non-dimensional 
sway  velocity  and  yaw  rate  are  all  set  to  zero  in  the  calculations.  This  avoids  division  by  zero 
errors  when  non-dimensionalizing  sway  velocity  and  yaw  rate. 

As  the  non-dimensional  forces  and  moments  are  multiplied  by  velocity  squared,  terms  less 
than  third  order,  e.g.,  r'v'  or  (v)2,  are  well  behaved  as  ship  speed  goes  to  zero.  The  coefficients 
for  third  order  and  higher  terms,  e.g.,  (r')2V,  are  linearly  reduced  to  zero  for  Fn  <  0.015. 

The  fraction  of  Tempest  rudder  forces  subtracted  from  the  hull  forces  is  linearly  reduced 
from  100%  to  0%  for  Fn  <  0.015.  While  this  appears  to  be  double  counting  the  rudder  at  near 
zero  speeds,  it  is  actually  gradually  switching  rudder  models  at  near  zero  speeds.  Tests  indicated 
the  coefficient  based  rudder  model  acted  as  though  there  was  no  rudder  at  near-zero  speeds. 

6.6  Rudder  and  Fin  forces 
6.6.1  Summary  of  Theory 

This  section  describes  the  theory  for  computing  the  rudder  force  for  a  ship  in  a  seaway. 

The  method  will  be  implemented  in  the  time  domain  dynamic  stability  code  Tempest.  The 
theory  is  based  primarily  on  the  empirical  fonnulae  developed  by  Wicker  and  Fehlner  (1958)  and 
is  developed  for  the  Level-0  implementation  of  Tempest. 

With  some  refinements,  the  theory  may  also  be  appropriate  for  the  higher  levels  of 
Tempest.  A  brief  synopsis  of  the  theory  will  be  described  in  this  section,  with  the  details 
provided  in  the  implementation  section  below.  The  ships  for  which  Tempest  is  currently 
envisioned  all  have  moveable  rudders,  as  opposed  to  flapped  rudders,  semi-balanced  rudders  on 
a  horn,  etc.  Therefore  the  rudder  model  described  here  applies  only  to  all  movable  rudders  and 
fins.  Prior  to  the  start  of  the  time  stepping  calculations  the  program  will  establish  fonnulae  for 
the  lift  and  drag  coefficients  of  the  rudder  as  a  function  of  angle  of  attack  based  on  empirical 
data.  Then  at  each  time  step  the  local  flow  velocity  at  the  rudder  will  be  computed.  The  velocity 
will  include  the  inflow  resulting  from  the  motion  of  the  ship,  the  wave  orbital  velocity,  and  the 
propeller  slipstream.  The  rudder  deflection  angle  will  be  provided  at  each  time  step  from  the 
autopilot  algorithm,  which  is  described  in  a  separate  document.  From  the  local  velocity  vector 
and  the  rudder  deflection  angle,  the  angle  of  attack  of  the  rudder  will  be  determined,  which  will 
be  used  to  compute  the  lift  and  drag  on  the  rudder.  The  lift  and  drag  forces  will  be  transformed 
to  forces  and  moments  acting  on  the  ship  in  the  ship-fixed  reference  frame.  If  the  rudder  is  only 
partially  submerged,  the  forces  will  be  adjusted  by  the  ratio  of  the  submerged  area  of  the  rudder 
to  the  total  area  of  the  rudder. 
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Figure  6-1  Definition  of  some  rudder  parameters.  (7  axis  into  page  on  left  half  figure,  X  axis  into 
page  on  right  half  figure.  The  value  of  </>r  is  negative  as  shown) 


6.6.2  Input 

The  user  will  be  required  to  provide  the  following  input  for  the  rudder.  The  first  set  of 
input  is  required  to  establish  the  position  and  geometry  of  the  rudder  and  the  second  set  of  input 
defines  some  of  the  coefficients  required  for  computing  the  rudder  force. 

The  position  and  geometry  of  the  rudder  will  be  computed  from  user  specified  values  for 
the  chord  length  of  the  rudder  at  the  root  and  tip,  and  the  (x,y,z)  offsets  of  the  quarter-chord  point 
at  the  root  and  tip  of  the  rudder:  c>,  ch  (xnynzr)  and  (xt,yt,z,)  shown  in  Figure  6-1. 

The  other  geometric  rudder  parameters  such  as  area,  span,  dihedral  angle,  etc.  can  all  be 
computed  from  these  values.  Internally  the  offsets  of  the  quarter-chord  point  at  the  tip  (xt,yt,zt) 
and  at  the  root  (xr,yr,zr)  will  be  referenced  to  the  ship-fixed  origin.  The  user  will  enter  the 
distances  to  these  points  forward  from  the  AP  and  above  the  baseline.  This  will  be  converted 
internally  to  the  position  of  the  quarter-chord  points  relative  to  the  ship  CG  in  the  ship-fixed 
frame. 

The  user  must  specify  several  other  parameters  to  compute  the  lift  and  drag  on  the  rudder. 
Some  values  will  have  preset  default  values,  in  which  case  the  user  will  have  the  option  to 
specify  a  value.  The  additional  inputs  are: 

•  ae-  the  effective  aspect  ratio  of  the  rudder.  The  geometric  aspect  ratio,  a,  of  the  rudder 
can  be  computed  as:  b /c  .  The  effective  aspect  ratio  accounts  for  the  mirror  image  effect 
when  the  rudder  is  flush  with  the  hull.  Typically  the  effective  aspect  ratio  is  assumed  to 
be  twice  the  geometric  aspect  ratio.  Flowever,  a  different  value  may  be  specified  if  a 
large  rudder  stool  is  used  or  the  rudder  has  a  large  tip  plate.  A  large  rudder  stool 
separates  the  rudder  from  the  hull  which  reduces  ae,  while  a  tip  plate  would  increase  the 
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loading  near  the  rudder  tip,  which  increases  ae .  The  effective  aspect  ratio  is  used  when 
computing  the  lift  and  drag  coefficients  and  the  stall  angle  of  the  rudder. 

•  Yr  -  flow  straightening  coefficient.  This  is  used  to  adjust  the  angle  of  attack  of  the  rudder 
to  account  for  the  flow  straightening  influence  of  the  hull  and  propeller 

•  Cdo  ~  minimum  sectional  drag  coefficient.  The  default  value  will  be  0.0065,  which 
corresponds  to  aNACA  0015  section 

•  Cdc  ~  crossflow  drag  coefficient  from  Whicker  and  Fehlner  (1958).  The  user  will  have 
the  option  to  either  enter  this  coefficient  directly,  or  to  specify  whether  the  tip  shape  is 
square  or  faired  and  have  the  program  compute  the  coefficient  based  on  the  data  in 
Whicker  and  Fehlner  (1958)  (which  expressed  the  coefficient  as  a  function  of  the  tip 
shape  and  taper  ratio). 

•  8—  rudder  deflection  angle.  This  will  be  provided  by  the  controller  algorithm  described 
in  Section  5  at  each  time  step. 

•  as  -  stall  angle.  The  program  will  compute  the  default  value  of  the  stall  angle  using  the 
formula  listed  in  “Step  1”  of  the  implementation  section.  The  user  may  choose  to 
override  the  default  value  by  specifying  a  value  directly 

•  Cat  -  Normal  force  coefficient  for  stalled  wing.  The  default  value  will  be  1.8  as 
recommended  by  Hoerner  (1975). 


6.6.3  Implementation 

Step  1  is  performed  before  the  start  of  the  time  stepping  procedure.  The  remaining  steps 
are  perfonned  at  each  time  step.  All  the  steps  will  be  perfonned  separately  for  each  rudder  on 
the  ship. 

Step  1.  Read  input,  compute  time  independent  values 

Prior  to  the  start  of  the  time  stepping  calculations  the  program  will  first  read  in  all  the  user 
supplied  inputs.  Then  the  geometry  information  required  for  computing  the  force  and  center  of 
pressure  for  the  rudder  will  be  computed. 

The  mean  chord  and  span  are  computed  as: 

c=^(ct+cr)  ,  b  -  ^(y, -yr)  2  +(zt -zrf  (6.73) 


Then  the  area  and  geometric  aspect  ratio  of  the  rudder  can  be  computed: 


Ar  —be 


=  b. 


(6.74) 

The  dihedral  angle  of  the  rudder  and  the  sweep  angle  of  the  quarter  chord  line  will  be 
computed  as: 


(f)R  -  arctan 


(yr~yt) 


A  =  arctan 


O,  ~xt) 


(6.75) 


The  spanwise  location  of  the  center  of  pressure  will  be  computed  by  assuming  an  elliptical 
spanwise  loading  distribution,  and  in  the  chordwise  direction  it  will  be  assumed  that  the  center  of 
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pressure  lies  on  the  quarter  chord  line.  With  these  assumptions  the  location  of  the  center  of 
pressure  is: 
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The  velocity  will  be  computed  at  each  time  step  at  the  mid-span  point  along  the  quarter 
chord  line: 
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The  velocity  will  first  be  computed  in  the  ship-fixed  reference  frame,  and  then  be 
transformed  to  a  frame  with  the  z  axis  parallel  to  the  rudder  stock,  as  indicated  by  the  axis 
system  marked  (x  \y  \z  ’)  in  Figure  6- 1 .  The  transformation  matrix  used  to  compute  velocities  in 
the  rudder  frame  from  velocities  defined  in  the  ship-fixed  frame  is: 
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where  (u,v,w)  is  the  velocity  vector  in  the  ship-fixed  reference  frame  and  (u  ’,v’,w  ’)  is  the 
velocity  in  the  rudder  frame.  The  transfonnation  matrix  is  orthogonal,  so  its  inverse  is  equal  to 
its  transpose. 


dC  / 

The  lift  slope  coefficient,  c  L/^a  ■>  and  the  stall  angle,  as,  can  be  computed  based  on  the 

geometry  and  user  specified  value  of  the  effective  aspect  ratio.  The  formula  used  for  the  lift 
slope  coefficient  is  that  developed  by  Whicker  and  Fehlner  (1958)  based  on  their  experimental 
data: 
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v  J 

If  a  value  for  the  cross  flow  drag  is  not  specified  by  the  user,  a  value  should  be  computed 
based  on  the  Whicker  and  Fehlner  (1958)  data: 


X  =  c/ 

A, 

CDC  -  0. 1  +  0.7/1  ,  for  rudders  with  faired  tips  (6.80) 

CDC  =  0. 1  + 1 .6/1  ,  for  rudders  with  square  tips 


Whicker  and  Fehlner  (1958)  recorded  the  stall  angles  for  the  fins  examined  in  their  wind 
tunnel  experiments  and  Lloyd  (1989)  fit  formulas  to  their  data.  The  rudder  will  operate  in  the 
slipstream  of  the  propeller,  and  the  turbulence  in  this  flow  delays  the  stall  to  a  higher  angle  from 
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what  was  recorded  in  the  Whicker  and  Fehlner  (1958)  wind  tunnel  tests.  If  a  value  for  stall  angle 
is  not  specified  by  the  user,  Tempest  will  use  the  formula  derived  by  Lloyd(1989)  to  specify  the 
rudder  stall  angle  with  an  additional  10  degrees  added  to  account  for  the  turbulence  in  the 
propeller  slipstream: 


as  =1.225 -0.445oe  +  0.075fi£  radians  forae<3.0 

as  =  0.565  radians  forae>3.0  (6.81) 


Step  2.  Compute  velocity  vector  at  midspan  point  on  the  rudder. 

At  each  time  step  the  velocity  vector  at  the  mid-span  point  along  the  rudder  quarter-chord 
line  (xm,ym,zm)  is  computed.  The  velocity  is  computed  by  first  adding  the  velocity  from  the 
motion  of  the  ship  to  the  wave  orbital  velocity,  and  then  adding  the  velocity  added  by  the 
propeller  slipstream.  The  procedure  for  computing  the  wave  orbital  velocity  at  an  arbitrary 
point  in  space  is  described  in  a  separate  document.  The  computations  for  the  wave  orbital 
velocity  are  performed  first  in  the  earth-fixed  frame,  but  must  be  transformed  to  the  ship-fixed 
frame  before  being  added  to  the  velocities  from  the  motion  of  the  ship.  The  velocity  vector  at 
the  rudder  mid-span  point  (xm,ym,zm)  is  computed  as: 
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(6.82) 


where:  ( u,v,w )  is  the  velocity  of  the  ship  defined  at  the  ship-fixed  origin,  (p,q,r )  is  the  rotational 
velocity  of  the  ship,  (uw,  vw,ww)  is  the  wave  orbital  velocity  at  the  center  of  the  propeller  disk,  vvF 
is  the  wake  fraction  coefficient,  which  is  discussed  in  Section  6.7,  and  UPROp  is  the  added  axial 
velocity  from  the  propeller.  All  the  above  quantities  are  defined  in  the  ship-fixed  reference 

frame.  The  negative  sign  before  the  first  and  last  vectors  in  the  formula  for  UR  indicates  that  the 
inflow  velocity  into  the  rudder  is  in  the  opposite  direction  from  the  ship  motion  and  the 
additional  velocity  from  the  propeller  slipstream  acts  along  the  negative  x  axis  in  the  ship  frame. 

The  additional  velocity  in  the  ship-fixed  x  direction  from  the  propeller  slipstream  is 
computed  following  the  method  described  in  Soding  (1981  and  1999).  First  the  axial  propeller 
slipstream  velocity  far  downstream  from  the  propeller  is  computed  from  the  momentum  equation 
as: 


^oo  —  U A +  CT  (6.83) 

From  the  propeller  force  calculations,  the  values  for  KT  and  J  are  available,  and  from  these 
CT  and  UA  can  be  computed  as: 
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(6.84) 


63 


The  radius  of  the  propeller  slipstream  far  downstream  from  the  propeller  based  on 
continuity  is: 


R,  =  R, 


V„+U, 


2F 


(6.85) 


where  R0  is  the  propeller  radius  equal  to  'AD.  The  radius  of  the  propeller  slipstream  at  the 
location  of  the  rudder  is  approximated  by  the  fonnula: 


Rx 
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(6.86) 


where  c  is  the  distance  from  the  propeller  to  the  rudder  normalized  by  the  propeller  radius: 
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The  velocity  in  the  propeller  slipstream  at  the  location  of  the  rudder  is  then  computed  from 
the  continuity  equation: 
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The  additional  velocity  from  the  propeller  is  added  to  the  rudder  inflow  velocity,  which  can  then 
be  computed  by  subtracting  the  inflow  velocity  to  the  propeller: 

U prop  =  VxRUD  ~ UA  (6.89) 

For  the  inclusion  of  the  propeller  slipstream  velocity  on  the  rudder,  the  inclination  of  the 
propeller  shaft  is  ignored  (at  least  for  Level-0).  Soding  (1999)  and  Brix  (1993)  also  include 
some  additional  corrections  to  the  influence  of  the  propeller  on  the  rudder  lift  by  accounting  for 
the  influence  of  turbulent  mixing  in  the  propeller  slipstream  and  the  influence  of  the  finite  width 
of  the  slipstream.  These  terms  will  not  be  included  in  the  Level-0  implementation  of  Tempest. 


Step  3.  Compute  angle  of  attack  for  the  rudder. 

To  compute  the  angle  of  attack  on  the  rudder,  the  inflow  velocity 
transfonned  into  the  rudder  frame  (the  (x  \y  ’,z  ’)  frame  shown  in  Figure 
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into  the  rudder  is  first 
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The  angle  of  attack  is  then  determined  by  examining  the  velocity  components  in  the  x  V 
plane  as  shown  in  Figure  6-2.  The  angle  of  the  inflow  into  the  rudder  relative  to  the  x  ’  axis  is 
computed  as: 


[3  =  arctan 


(6.91) 


This  angle  is  adjusted  by  the  flow  straightening  coefficient,  yR,  to  account  for  the  flow 
straightening  influence  of  the  hull  and  propeller.  Molland  and  Turnock  (2002)  discussed  the  use 
of  the  flow  straightening  coefficient  and  performed  a  set  of  experiments  to  determine  values  of  yR 
for  several  simple  configurations.  In  Tempest,  the  user  will  provide  a  value  of  yR  as  input,  with 
the  default  value  being  1.0. 

P'  =  YrP  (6  92) 

The  angle  of  attack  is  then  computed  by  subtracting  /?’  from  the  rudder  deflection  angle,  c\ 
which  is  provided  by  the  autopilot  routine  at  each  time  step: 


a  =  ksS  -  P' 


(6.93) 


The  rudder  deflection  angle  fiis  defined  as  positive  using  a  right  hand  rule  about  the  z'  axis, 
which  points  upwards  (out  of  the  paper  in  Figure  6-2)  resulting  in  a  turn  to  starboard  as  shown  in 
Figure  6-2.  The  positive  direction  for  the  angles  a  and  P are  defined  in  the  same  manner.  The  Ay- 
factor  accounts  for  fixed  portions  of  the  rudder  (the  rudder  stool)  that  effectively  reduce  the 
angle  of  attack.  k§  is  set  equal  to  the  ratio  of  the  planform  area  of  the  movable  rudder  over  the 
total  planfonn  area  of  the  rudder  and  stool.  For  small  to  moderate  rudder  deflection  angles,  the 
rudder  and  stool  are  modeled  as  a  single  lifting  surface  with  a  variation  in  angle  of  attack  along 
the  span. 


Figure  6-2  Inflow  into  the  rudder  and  calculation  of  angle  of  attack. 
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Step  4.  Compute  the  Lift  and  Drag  on  the  rudder 

For  cases  where  the  angle  of  attack  is  below  the  stall  angle,  the  lift  and  drag  coefficients  are 
computed  from  the  empirical  fonnulae  derived  in  Whicker  and  Fehlner  (1958). 

For  \a\  <  as  : 
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where - and  Cdc  were  computed  in  Step  1.  For  cases  where  the  angle  of  attack  is  greater 

da 

than  the  stall  angle,  it  will  be  assumed  that  the  flow  is  completely  separated  at  the  leading  and 
trailing  edge  of  the  rudder,  and  the  force  is  normal  to  the  chord  line  of  the  rudder.  The  resultant 
lift  and  drag  coefficients  in  this  condition  are  computed  from  the  formulae  found  in  Hoerner 
(1975): 

For  led  >as  , 


Cl=B-Cn  sinacosa  +  C,  max  ■  BQ 


(6.95) 
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where  B,  B0,  Cl, max,  and  Co.max  have  been  included  to  create  a  cubic  transition  between  the  fully 
stalled  flat  plate  and  the  maximum  Cl  condition  when  \a\  >  as  and  are  defined  as: 


for  \a\>\25-as  : 
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After  the  lift  and  drag  coefficients  are  computed,  the  lift  and  drag  can  be  calculated.  The 
velocity  used  to  compute  the  forces  from  the  lift  and  drag  coefficients  is  the  component  of  the 
inflow  velocity  to  the  rudder  in  the  x  ’y  ’  plane. 
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(6.98) 


The  lower  case  d  is  used  here  for  drag  to  avoid  confusion  with  the  symbol  used  earlier  for 
the  propeller  diameter. 


Step  5.  Correct  rudder  lift  and  drag  for  partial  rudder  submergence. 

If  the  rudder  is  not  completely  submerged,  the  lift  and  drag  on  the  rudder  will  be  corrected 
to  account  for  this.  The  lift  and  drag  will  be  adjusted  based  on  the  approximate  percentage  of  the 
rudder  area  that  is  submerged.  This  is  equivalent  to  replacing  AR  in  the  fonnulae  for  lift  and  drag 
listed  above  with  the  submerged  rudder  area  A  RS.  To  approximate  A  KS  the  wave  elevation  at  the 
rudder  is  first  computed.  Let  zWR  be  the  wave  elevation  above  the  rudder  location  (xm,ym).  If  the 
rudder  tip  is  above  the  wave:  z,  >  zWR,  then  the  rudder  is  completely  out  of  the  water  and  the 
rudder  lift  and  drag  are  set  to  zero.  If  the  rudder  root  is  below  the  wave:  zr  >  zWR,  then  the  rudder 
is  completely  submerged  and  no  adjustment  to  the  lift  and  drag  from  step  4  is  required. 

Otherwise  the  submerged  area  of  the  rudder  is  computed  as  follows.  First  compute/  the 
percentage  of  the  rudder  span  that  is  submerged: 
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The  submerged  area  of  the  rudder  can  then  be  approximated  as: 
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(6.100) 


The  lift  and  drag  on  the  rudder  are  then  recalculated  using  ARS  in  place  of  AR. 

L  ~  P  ^Rs^L  (6.101) 

d  —  ^  P  V r  ^rs  C d 


Step  6.  Compute  rudder  forces  and  moments  in  the  ship-fixed  frame. 

The  lift  and  drag  force  on  the  rudder  act  in  the  x  ’y  ’  plane  in  the  direction  perpendicular  and 
parallel  to  the  inflow  velocity  to  the  rudder  defined  by  the  angle  J3’.  To  resolve  the  forces  in  the 
ship-fixed  frame,  the  forces  are  first  transformed  to  forces  in  the  x  ’  and  y  ’  directions,  and  then 
are  transformed  to  the  ship-fixed  frame.  The  forces  in  the  x  ’  and y  ’  directions  are: 

Fxr  =~d  cos  P'-L  sin  /?'  (6.102) 

Fyr  =  L  cos  P'-d  sin  /?' 
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The  forces  are  then  transformed  to  the  ship-fixed  frame  using  the  transpose  of  the 
transformation  matrix  defined  in  Step  1 . 
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(6.103) 


The  rudder  force  acts  at  the  center  of  pressure,  ( xc,yc,zc ),  which  was  computed  in  Step  1 . 
The  moments  which  result  from  the  rudder  force  are  computed  from  FR  and  ( xc,yc,zc ): 
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Fr  and  M R  are  the  force  and  moment  from  the  rudder  acting  at  the  ship-fixed  reference 
frame  origin,  which  is  the  ship  CG  in  Level-0.  These  two  vectors  are  added  to  the  other  force 
components  to  compute  the  total  force  and  moment  action  on  the  ship  at  each  time  step. 


6.6.4  Notes  on  ship  specific  empirical  data 

Rudder  data  either  from  measurements  in  the  Large  Cavitation  Tunnel  (LCC)  or  from  CFD 
computations  are  available  for  the  rudder  designs  for  some  ships.  When  available  this  data 
should  be  used  to  help  specify  the  lift  and  drag  coefficients  more  accurately  for  the  rudder  model 
implemented  in  Tempest.  Within  the  Level-0  theory  discussed  in  this  document,  the  data  from 
LCC  tests  could  be  used  to  specify  the  stall  angle  and  to  determine  a  value  for  the  effective 
aspect  ratio  that  closely  matches  the  lift  slope  curve. 

One  difficulty  in  basing  the  lift  and  drag  coefficients  solely  on  the  data  from  the  LCC  tests, 
is  that  these  tests  have  been  performed  only  with  the  ship  traveling  straight  ahead  with  the  rudder 
deflected  at  different  angles.  The  influence  of  the  ship  drift  angle  and  cross  flow  from  the  wave 
induced  velocity  would  need  to  be  accounted  for  in  some  manner.  If  the  stool  on  which  the 
rudder  sits  has  a  longer  span  than  a  typical  rudder  stool,  the  lift  on  the  rudder  stool  should  be 
considered.  The  rudder  force  model  currently  computes  the  force  only  on  the  rudder,  not  on  the 
stool.  The  Level-0  maneuvering  force  model  was  developed  from  rotating  ann  test  data  with  the 
rudders  and  rudder  stool  in  place,  and  the  rotating  arm  tests  measured  the  total  force  on  the 
model  including  the  forces  from  rudder  and  the  rudder  stool.  The  force  on  the  rudder  is  then 
computed  and  subtracted  out  to  obtain  force  coefficients  for  building  the  maneuvering  model. 
However  the  force  from  the  rudder  stool  should  still  be  included  in  the  maneuvering  force  model, 
so  no  separate  calculation  of  the  force  from  the  rudder  stool  should  be  required,  if  the  Level-0 
maneuvering  model  is  used. 

The  stall  angle,  as,  can  be  estimated  from  empirical  data.  The  stall  angle  is  dependent 
primarily  on  the  aspect  ratio,  profile  shape,  thickness,  Reynolds  number,  surface  roughness,  and 
the  turbulence  in  the  inflow  (Brix  1993).  This  makes  an  accurate  prediction  of  the  stall  angle 
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difficult.  Figure  6-3  shows  the  empirical  formula  used  in  Level-0  Tempest  which  is  based  on 
Lloyd’s  analysis  of  the  Whicker  and  Fehlner  data. 


Stall  Angle  as  function  of  Aspect  Ratio 


Figure  6-3  Stall  angle  formula  derived  by  Lloyd  (1989)  based  on  Whicker  and  Fehlner  (1958) 


Several  empirical  formulae  are  available  for  computing  the  lift  coefficient  slope  ( 
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as  a  function  of  effective  aspect  ratio.  Figure  6-4  compares  the  formulae  proposed  by  Hoerner 
(1975),  Soding  (1999)  and  Whicker  and  Fehlner  (1958).  There  is  much  less  empirical  data 
available  for  the  lift  and  drag  on  rudders  at  angles  of  attack  significantly  over  the  stall  angle. 
Figure  6-5  is  taken  from  Hoerner  (1975)  and  show  the  lift  from  several  2-D  foil  experiments. 
Figure  6-6  shows  the  computed  lift  and  drag  coefficients  using  the  formulae  specified  for  the 
Tempest  Level-0  theory  over  the  full  range  of  angles  of  attack  from  0°  to  180°  for  a  rectangular 
foil  with  an  aspect  ratio  of  2.  The  curve  for  the  lift  coefficient  compares  well  with  the  2-D  test 
data  from  Hoerner  (1975). 
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Lift  Slope  as  a  function  of  Aspect  Ratio 


Figure  6-4  Comparison  of  formulae  for  lift  coefficient  slope. 


figure  33.  Variation  of  lift  coefficient  for  angles  of  attack  0  to 
1 80  degrees. 


Figure  6-5  Lift  on  2D  foils,  past  stall  angle,  reproduced  from  Hoemer  (1975) 
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Figure  6-6  Lift  and  Drag  coefficients  computed  for  a  rectangular  foil  with  an  aspect  ratio  of  2.0, 
based  on  the  formulae  described  in  the  Tempest  Level-0  theory. 


6.7  Propulsor  Force 
6.7.1  Summary  of  Theory 

This  section  describes  the  theory  for  computing  the  propeller  force  for  a  ship  in  a  seaway. 
This  method  will  be  implemented  in  the  time  domain  dynamic  stability  code  Tempest.  The 
theory  is  initially  developed  for  the  Level-0  implementation  of  Tempest. 

With  some  refinements,  the  theory  may  also  be  appropriate  for  the  higher  levels  of 
Tempest.  The  theory  is  developed  initially  for  conventional  fixed-pitch  propellers.  Extensions 
to  the  theory  would  be  required  to  model  controllable  pitch  propellers,  waterjets,  and  other  non- 
conventional  propulsors. 

A  brief  synopsis  of  the  theory  will  be  described  in  this  section,  with  the  details  provided  in 
the  implementation  section  below.  The  propeller  thrust  will  be  determined  by  interpolating  the 
user-supplied  open  water  propeller  thrust  curve  (KT  vs.  J)  at  each  time  step.  The  advance 
coefficient,  J,  will  be  computed  at  each  step  from  the  instantaneous  motion  of  the  ship  and  the 
wave  orbital  velocity.  The  velocity  will  be  computed  at  a  single  point  located  at  the  center  of  the 
propeller  disk.  The  velocity  vector  at  this  point  will  be  decomposed  into  a  component  parallel  to 
the  propeller  shaft  and  a  cross  flow  velocity  component.  The  influence  of  the  ship  hull  on  the 
propeller  inflow  will  be  accounted  for  using  the  Taylor  wake  fraction  coefficient  and  the  increase 
in  hull  resistance  from  the  propeller  will  be  accounted  for  using  a  thrust  deduction  fraction. 
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Propeller  forces  that  are  perpendicular  to  the  shaft  axis  resulting  from  cross  flow  will  be 
estimated  using  empirical  formulae  described  by  Fuhs  and  Dai  (2007).  If  the  propeller  is  only 
partially  submerged,  the  forces  will  be  adjusted  by  the  ratio  of  the  submerged  area  of  the 
propeller  disk  to  the  total  area  of  the  disk. 


6.7.2  Input 

The  user  will  be  required  to  provide  the  following  input  to  characterize  the  propeller 
performance.  The  user  will  supply  a  set  of  points  describing  the  open  water  propeller  thrust 
curve:  KT  vs.  J.  The  coefficients  KT  and  J  are  defined  below: 
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The  curve  should  be  defined  by  a  sufficient  number  of  points  for  KT  to  be  interpolated  at 
each  time  step,  and  should  include  the  entire  first  quadrant  of  the  open  water  curve  including  a 
value  for  J  =  0  up  to  the  value  of  J  where  KT  =  0.  Negative  values  of  J  may  also  be  entered,  to 
allow  for  reverse  flow  into  the  propeller. 

In  addition  the  user  will  supply  the  following  values: 

D  diameter  of  propeller 

rpm  rotations  per  minute  (converted  to  rotations  per  second  internally) 
p  density  of  water 

wp  Taylor  wake  fraction  coefficient 

tp  thrust  deduction  coefficient 

xp,yp,zp  position  of  the  center  of  the  propeller  disk  in  ship-fixed  frame 

(The  user  will  enter  the  distance  of  the  propeller  forward  from  the  AP  and  above 
the  baseline,  this  will  be  converted  internally  to  the  position  of  the  propeller 
relative  to  the  ship  CG  in  the  ship-fixed  frame.) 

as,  Ps  angles  defining  propeller  shaft  orientation 

y  flow  straightening  factor 

RHP  integer  flag  to  indicate  whether  propeller  is  right  hand  turning  (+1)  or  left  hand 
turning  (-1).  A  right  hand  turning  propeller  will  turn  clockwise  when  viewed 
from  the  stern  looking  forward. 


The  Taylor  wake  fraction  coefficient  is  used  to  adjust  the  velocity  into  the  propeller  to 
account  from  the  presence  of  the  ship  hull  and  is  defined  as: 
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where  V  is  the  axial  velocity  at  the  location  of  the  propeller  in  the  absence  of  the  ship  hull  and  Va 
is  the  axial  component  of  local  velocity  vector  at  the  center  of  the  propeller  disk  in  the  presence 
of  the  hull. 

The  thrust  deduction  fraction  describes  the  additional  resistance  on  the  hull  resulting  from 
the  propeller  -  hull  interactions.  Traditionally  this  additional  resistance  is  thought  of  as  a 
reduction  of  the  propeller  thrust.  The  thrust  deduction  fraction  is  defined  as: 

t  T~Rhun  (6.108) 

T 

where  Rhuii  is  the  towed  resistance  of  the  hull  without  the  propeller  and  T  is  the  propeller  thrust. 
Methods  for  determining  wp  and  tp  for  twin  screw  destroyers  can  be  found  in  Hurwitz  (1980). 

The  flow  straightening  factor  is  an  optional  input  variable  that  will  be  used  to  adjust  the 
magnitude  of  the  cross  flow  velocity  to  account  for  the  influence  the  ship  hull  and  skeg  have  on 
straightening  the  flow  into  the  propeller.  It  is  defined  as: 

(6.109) 

where  </>s  is  the  angle  between  the  total  velocity  vector  at  the  propeller  and  the  shaft  axis  in  the 
absence  of  the  hull  and  (f>s  is  the  same  angle  accounting  for  the  presence  of  the  hull.  A  default 
value  of  7=  1.0  will  be  used  if  the  user  does  not  supply  a  value. 


6.7.3  Implementation 

Step  1  is  performed  before  the  start  of  the  time  stepping  procedure.  The  remaining  steps 
are  performed  at  each  time  step.  All  of  the  steps  will  be  performed  separately  for  each  propeller 
on  the  ship. 


Step  1.  Read  input,  compute  time  independent  values 

Prior  to  the  start  of  the  time  stepping  calculations  the  program  will  first  read  in  all  the  user 
supplied  inputs.  A  unit  vector  defining  the  direction  of  the  propeller  shaft  in  the  ship-fixed 
reference  frame  will  be  defined  using  the  fonnula: 


n 


s 


^ -  cos  J3S  COS  <2^ 

sin  Ps 

-  sin  as  cos  J3S 


(6.110) 


where,  as  describes  the  vertical  inclination  of  the  propeller  shaft  and  is  defined  as  the  angle 
between  the  x-axis  of  the  ship-fixed  frame  and  the  projection  of  the  shaft  axis  onto  the  vertical 
plane  of  symmetry  of  the  ship.  A  positive  angle  indicates  that  the  shaft  points  downwards.  J3s 
describes  the  orientation  of  the  propeller  shaft  in  the  horizontal  plane  relative  to  the  ship  x-axis 
with  a  positive  angle  indicating  that  the  shaft  points  to  port. 
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Step  2.  Compute  velocity  vector  at  center  of  propeller  disk 


At  each  time  step  the  velocity  vector  at  the  center  of  the  propeller  disk  is  computed.  The 
velocity  is  computed  by  adding  the  velocity  from  the  motion  of  the  ship  to  the  wave  orbital 
velocity.  The  procedure  for  computing  the  wave  orbital  velocity  at  an  arbitrary  point  in  space  is 
described  in  a  separate  document.  The  computations  for  the  wave  orbital  velocity  are  perfonned 
first  in  the  earth-fixed  frame,  but  must  be  transformed  to  the  ship-fixed  frame  before  being  added 
to  the  velocities  from  the  motion  of  the  ship.  The  velocity  vector  at  the  center  of  the  propeller 
disk  (xp,yp,zp)  is  computed  as: 


"  u  +  zpq-ypr  ' 

fu  } 

V  +  Xpr~ZpP 

+ 

yw+ y  pP  ~  x  pq  j 

(6.111) 


where,  ( u,v,w )  is  the  velocity  of  the  ship  defined  at  the  ship-fixed  origin,  (p,q,r )  is  the  rotational 
velocity  of  the  ship,  (xp,yp,zp)  is  the  position  of  the  center  of  the  propeller  disk  and  (uw,vw,ww)  is 
the  wave  orbital  velocity  at  the  center  of  the  propeller  disk.  All  the  above  quantities  are  defined 

in  the  ship-fixed  reference  frame.  The  negative  sign  before  the  first  vector  in  the  formula  for  U p 
indicates  the  inflow  velocity  into  the  propeller  is  in  the  opposite  direction  from  the  ship  motion. 


Step  3.  Decompose  velocities  into  components  parallel  and  perpendicular  to  shaft. 
Compute  cross  flow  angle  and  velocity. 

The  component  of  the  velocity  parallel  to  the  propeller  shaft  is  the  dot  product  of  the 
velocity  vector  with  the  unit  vector  parallel  to  the  shaft  axis: 

U ,,  =  UP-ns  =cosfis  cosas(u  +  zpq-ypr)-smps(v+xpr-zpP)  +  sinas  cos fis(w+ypp-xpq) 


(6.112) 

The  cross  flow  velocity  is  then  the  difference  between  the  total  velocity  vector  and  the 
component  parallel  to  the  shaft: 

U  cross  =  Up  ~  U  Ans  (6. 113) 

The  unit  normal  vector  in  the  direction  of  the  cross  flow  can  be  computed  from  Ucross : 


U, 


CROSS 


U, 

and  the  cross  flow  angle,  fa,  is  defined  as: 

<j>s  =  arctan 


CROSS 


f  TT 

/  \ 

U  CROSS 

/ 

VUA 

V 

y 

(6.114) 


(6.115) 


74 


0s  should  always  be  a  positive  angle  between  0  and  n/2.  The  unit  vector  defining  the  direction 
perpendicular  to  both  the  shaft  axis  and  the  cross  flow  is  computed  as  the  cross  product  of  the 
unit  vectors  defining  the  shaft  axis  and  cross  flow  directions: 


nN  =nsxnc  = 


nS2nCi  nS3nC2 


\nS\nC2  nS2nC\  ) 


(6.116) 


Defined  in  this  way,  nN  will  point  upwards  when  the  cross  flow  is  to  starboard  and  point 
downwards  when  the  cross  flow  is  to  port.  The  magnitude  of  the  vector  formed  by  the  cross 


product  of  two  vectors  is  defined  as: 


a  xb 


=  \a\ 


sin/ ,  where  6  is  the  angle  between  the  two 


vectors.  In  the  case  of/7s  x  nc ,  both  vectors  are  unit  vectors  and  they  are  perpendicular,  so  the 

angle  between  them  should  be  90°.  They  should  also  have  a  magnitude  of  the  cross  product 
vector  equal  to  one.  Computing  the  magnitude  of  this  vector  would  provide  a  useful  check. 


Step  4.  Compute  the  propeller  thrust,  side  force  and  normal  force 

The  propeller  thrust  is  computed  by  interpolating  the  user  supplied  KT  vs.  /.  The  advance 
coefficient  /  is  first  computed  from  the  axial  velocity  adjusted  to  account  for  the  presence  of  the 
hull  using  the  Taylor  wake  fraction: 

j_UA(\-wF )  (6.117) 

nD 


The  value  for  Kj  is  then  obtained  by  interpolating  the  user  supplied  data.  If  the  input  data 
does  not  cover  this  range,  it  should  be  extrapolated  to  include  these  end  points,  with  a  warning 
printed  to  the  output  file.  Negative  values  of  J  may  also  be  specified  to  provide  data  for  reverse 
flow  into  the  propeller,  where  UA  is  negative  and  n  is  positive.  If  the  computed  value  for  /  is 
outside  the  range  of  user-supplied  data  for  KT  vs.  /,  extrapolation  should  not  be  used,  and  the 
value  used  for  Kj  should  remain  constant  outside  the  range  of  specified  data.  In  the  Level-0 
implementation  it  is  assumed  that  the  revolutions  per  second,  n,  is  a  constant  value  input  by  the 
user.  Later  implementations  may  improve  on  this  to  account  for  the  influence  of  the  propeller 
loading  on  the  rpm  of  the  propeller.  If  the  value  for  the  rotational  speed,  n,  is  zero,  then  both  / 
and  Kt  should  be  set  to  zero  along  with  all  the  propeller  forces.  The  KT  vs.  /  should  cover  values 
of  /ranging  from  0  to  the  value  of  /where  the  thrust  is  zero. 


After  computing  KT,  the  thrust  is  computed  as: 

T  =  Kt  p  n  D4 


(6.118) 


For  cases  where  the  UA  is  very  small,  the  cross  flow  angle,  0s,  as  computed  above  may  give 
values  close  to  90°  even  for  small  cross  flow  velocities.  It  would  be  expected  that  the  induced 
velocity  into  the  propeller  would  dominate  the  local  flow  in  these  cases.  Therefore,  after  the 
thrust  is  computed,  the  velocity  induced  by  the  propeller  at  the  center  of  the  propeller  disk  will 
be  computed  from  momentum  theory  as: 
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(6.119) 


U 


A  PROP 


If  the  value  computed  for  U a  prop  is  greater  than  the  value  of  UA,  then  the  cross  flow  angle, 
</>s,  will  be  recomputed  replacing  UA  with  UA  prop  in  the  formula  listed  in  Step  3. 

Next  the  side  force  and  normal  force  are  computed  based  on  the  empirical  formulae  derived 
in  Fuhs  and  Dai  (2007). 

Fcross  =  1.716  T  sin(/$s)  ^  120) 

Fnorm  =0.2  Fcross 


where  /is  an  optional  input  value  with  a  default  value  of  1.0  that  may  be  used  to  account  for  the 
flow  straightening  influence  of  the  hull. 


Step  5.  Compute  submerged  area  of  the  propeller  disk 

The  following  is  one  way  to  calculate  the  submerged  area  of  the  propeller  disk.  The 
calculation  of  the  submerged  area  of  the  propeller  disk,  As  ,  is  perfonned  in  the  earth-fixed  frame. 

The  vertical  position  of  the  center  of  the  propeller  disk  in  the  earth-fixed  frame,  Zp  ,  is  first 
computed  from  ( xp,yp,zp )  using  the  transformation  formulae  described  in  the  Section  3.2.  The 
“E”  superscripts  on  variables  in  this  section  indicate  that  the  values  are  defined  in  the  earth-fixed 

frame.  Zw  is  the  wave  elevation  directly  above  the  propeller.  It  is  assumed  that  the  wave  is 

long  relative  to  the  diameter  of  the  propeller,  so  the  wave  height  is  assumed  to  be  constant  above 
the  propeller.  The  inclination  of  the  shaft  is  also  ignored  during  the  calculation  of  the  submerged 
area.  The  total  area  of  the  propeller  disk  is:  AD  -  kR1  ,  where  R  is  the  propeller  radius,  Vi  D. 


If  the  value  of  ZEV  is  greater  than  Z'p  +  R ,  then  the  propeller  is  completely  submerged  and 

As  -  Ad.  If  the  value  of  is  less  than  ZEp  —R ,  then  the  propeller  is  completely  out  of  the 

water  and  A  s  =  0  .  For  all  other  conditions,  the  propeller  is  partially  submerged  and  the 

submerged  area  must  be  computed.  This  calculation  can  be  made  by  first  computing  the  angle  6 
shown  in  Figure  6-7. 


9  =  arcsin 


f  ry  E 

^ w 


'  E  \ 


R 


(6.121) 


The  angle  9  will  be  between  and  +  •  It  will  t>e  positive  if  the  wave  is  higher  than 

the  center  of  the  propeller  disk  and  negative  if  the  wave  is  below  the  center  of  the  propeller  disk. 
The  submerged  area  can  now  be  computed  in  terms  of  6  and  R  as: 

6  a  7 

cosf?sin0H - h  0 

2  , 


AS=R- 


(6.122) 
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Figure  6-7  Calculation  of  submerged  area  of  propeller  disk 

An  alternate  way  to  approximate  the  area  is  to  panelize  the  propeller  disk  and  trim  the 
panels  using  the  undisturbed  free  surface  as  a  cutting  surface.  This  is  the  approach  taken  by 
Tempest. 


Step  6.  Compute  forces  and  moments  in  ship-fixed  frame. 

The  thrust  is  in  the  direction  of  the  shaft  axis,  but  in  the  opposite  direction  from  the  unit 
vector  defining  the  shaft  axis  direction,  so  the  thrust  is  applied  in  the  direction:  -  ns  .  The  side 

force  F cross  is  in  the  direction  of  the  cross  flow  nc  for  both  right  hand  turning  and  left  hand 
turning  propellers.  The  direction  of  the  normal  force,  FNORM,  depends  on  both  the  direction  of 
the  cross  flow  and  whether  the  propeller  is  right  hand  turning  or  left  hand  turning.  For  a  right 
hand  turning  propeller,  the  normal  force  will  point  upwards  if  the  cross  flow  is  to  port  and 
downwards  if  the  cross  flow  is  to  starboard.  The  nonnal  force  will  point  in  the  opposite  direction 
for  a  left-handed  propeller.  The  unit  vector  nN  is  defined  such  that  it  points  upwards  when  the 

flow  is  to  starboard.  Therefore  the  direction  of  the  nonnal  force  can  be  defined  by  the  vector: 

-  RIIS  nc ,  where  RHS  is  + 1 .0  for  a  right  hand  propeller  and  - 1 .0  for  a  left  hand  propeller. 

The  propeller  force  vector  in  the  ship-fixed  reference  frame  can  be  computed  from  the 
force  components  parallel  and  perpendicular  to  the  propeller  shaft: 

F  =—(-Tn  +F  ft  -F  RHPn 

i  p  y  i  Us  t  i  CROSSnc  1  CROSS  "n 

ad 

The  x  component  of  the  force  in  the  ship-fixed  frame  is  adjusted  to  account  for  the 
additional  drag  on  the  hull  resulting  from  the  interaction  between  the  hull  and  propeller  using  the 


U  ~tp) 

1 

1 


(6.123) 
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user  specified  thrust  deduction  fraction,  tP.  All  three  components  of  the  force  are  multiplied  by 


the  ratio  of  the  submerged  disk  area  to  the  total  disk  area: 


The  moments  which  result  from  the  propeller  force  are  computed  from  Fp  and  (xp,yp,zp): 


FpZyP  f„Zi 


Fpxzp  Fpzxp 
yFpYXp  —  Fpxy  p  y 


(6.124) 


Fp  and  Mp  are  the  force  and  moment  from  the  propeller  acting  at  the  ship-fixed  reference 
frame  origin,  which  is  the  ship  CG  in  Level-0.  These  two  vectors  are  added  to  the  other  force 
components  to  compute  the  total  force  and  moment  action  on  the  ship  at  each  time  step. 
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